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Abstract 

In  the  first  part  of  this  work  the  overland  flow  production  mechanism  used  in  the 
Willgoose-Bras-Rodriguez-Iturbe  channel  network  and  catchment  evolution  model  was 
modified  from  the  Hortonian  to  the  subsurface  saturation  runoff  mechanism.  Two  important 
differences  were  found  oin  the  behavior  of  the  new  model;  one  was  in  the  evolution  of  the 
hypsometric  curves  of  the  catchment;  the  second  was  in  the  importance  of  mass  movement 
(e.g.  creep  and  landsliding)  sediment  transport  in  the  overall  behavior  of  the  system.  A  new 
nondimensional  number  related  to  the  threshold  of  saturation  in  hillslopes  was  added  to  the 
list  of  nondimensional  numbers  belonging  to  the  original  model.  Common 
geomorphological  statistics  were  measured  on  the  simulated  catchments  and  found  to  be 
similar  to  field  data. 

The  second  part  of  this  work  examines  the  sensitivity  to  initial  conditions  present  in  the 
WBR  model.  Using  an  appropiate  measure  between  elevation  fields,  the  evolution  of 
different  catchments  was  examined.  Exponential  separation  in  phase  space  between 
trajectories  of  the  system  was  found  in  a  variety  of  tests.  This  exponential  separation  is  the 
reason  for  the  apparent  randomness  present  in  die  evolution  of  the  model.  The  influence  of 
different  parameters  of  the  model  on  this  exponential  separation  was  also  examined. 
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Chapter  1 
Introduction 


1.1  Scope 

Hydrologists  have  always  been  interested  in  the  geomorphologic  aspects  of 
catchments  and  river  networks.  An  understanding  of  the  relationships  between  basin  shapes 
and  hydrologic  responses  is  fundamental  for  the  purpose  of  hydrologic  predictions, 
especially  in  ungaged  basins. 

Another  goal  is  to  try  to  understand  not  only  the  properties  and  the  response  of  a 
catchment  at  a  certain  rime  but  to  study  the  evolution  of  the  catchment  over  a  long  period  of 
time.  This  approach  has  potential  applications  in  the  study  of  geomorphologic  impacts 
produced  by  global  climate  changes.  Water,  as  a  vehicle  of  sediment  transport,  plays  a 
fundamental  role  in  the  process.  Recently  a  model  developed  by  Willgoose.  Bras  and 
Rodriguez- Iturbe  [47]  quantified  the  evolution  of  the  channel  network  and  the  catchment. 
This  work  uses  that  model  to  address  two  questions:  what  is  the  influence  of  the  overland 
flow  production  mechanism  in  the  overall  behavior  of  the  basin  and  why  is  there  an 

apparent  randomness  present  in  the  evolution  process  even  if  the  system  is  completely 
deterministic. 

12  Outline 

Chapter  2  presents  a  general  review  of  previous  work.  This  chapter  gives  basic 
definitions  and  describes  some  channel  network  models.  The  problems  described  in  the 
above  section  will  be  discussed  in  two  separate  chapters.  Chapter  3  studies  a  modification 
to  the  WBR  model  where  a  subsurface  saturation  mechanim  is  used  to  calculate  the 
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overland  flow  distribution.  Several  simulations  are  presented  and  the  geomorphological 
properties  of  the  catchment  and  networks  are  studied.  Also,  comparisons  between  the 
original  model  and  its  modified  version  are  made.  A  non-dimensional  form  of  the  model  is 
developed  which  is  useful  for  comparisons  among  basins  of  different  spatial  characteristics. 
The  influence  of  a  probabilistic  distribution  of  rainfall  over  the  overland  flow  distribution 
and  the  behavior  of  the  model  is  studied  in  the  last  section  of  Chapter  3. 

Chapter  4  addresses  the  hypothesis  of  the  existence  of  transient  chaos  in  the  system 
and  its  implication  in  the  system’s  sensitivity  to  perturbations  in  the  elevation  field,  a 
reflection  of  the  chaotic  nature  of  the  system.  The  possibility  of  the  existence  of  chaos 
would  explain  the  apparent  randomness  in  the  resulting  channel  networks. 

Chapter  5  summarizes  this  work  and  suggests  further  avenues  for  research.  Appendix 
A  presents  the  values  of  the  non-dimensional  parameters  of  the  different  simulations. 
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Chapter  2 

Literature  review 


2.1  Basic  Concepts 


The  flow  of  water  shapes  river  basins  but  at  the  same  time  it  is  the  form  of  the  basin 
that  determines  in  what  direction  water  flows.  These  two  complementary  effects  constitute 
the  link  between  hydrology  and  geomorphology:  hydrologic  behavior  of  basins  and. 
evolution  of  catchments  and  river  networks. 


The  quantitative  study  of  channel  networks  began  with  two  papers  by  Horton  [13], 
[14]  in  which  an  ordering  system  of  networks  was  introduced.  Later,  Strahler  [41]  modified 
this  system  and  it  is  his  ordering  that  is  most  commonly  used  today.  In  this  scheme  a 
stream  that  has  a  source  as  one  of  its  extreme  points  is  assigned  order  1.  Where  two  streams 
of  order  i  join,  the  channel  downstream  is  of  order  i+1.  If  two  channels  of  unequal  order 
join,  the  channel  downstream  has  the  same  order  as  the  higher  order  incoming  stream.  The 
order  of  the  outlet  stream  is  the  basin  order,  £2.  Exterior  links  arc  those  between  a  source 
and  a  junction  and  interior  links  are  those  between  two  junctions.  Various  link  properties 
like  length,  area  and  slope  are  of  interest  to  hydrologists.  One  of  the  purposes  of  network 
ordering  is  to  find  general  relationships  between  links  of  different  orders  across  the  basin. 


The  most  common  geomorphological  measures  are  the  Horton/Strahler  statistics 
which  are  stream  based  are: 


(2.1) 


(22) 
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=  (2-3) 

^oi-l 


where  Nu  is  the  number  of  streams  of  order  co,  is  the  mean  length  of  streams  of  order  co, 
is  the  mean  area  contributing  to  streams  of  order  co,  and  Sw  is  the  mean  slope  of  streams 
of  order  co.  Rfe,  R;,  Ra  and  R^  are  the  bifurcation,  length,  area  and  slope  Horton's  ratios.  The 
interest  of  these  ratios  is  that  they  are  approximately  constant  when  plotted  against  order,  a 
property  that  is  usually  called  Horton’s  law.  A  lot  of  work  has  been  devoted  to  the  analysis 
of  these  ratios  in  real  networks  under  a  variety  of  conditions.  These  catchment  statistics  can 
help  to  understand  hydrologic  properties  of  the  basin  and  the  work  on  the  geomorphologic 
instantaneous  unit  hydograph  GIUH  by  Rodriguez-Iturbe  and  Valdes  [31]  and  Gupta  et  al. 
[10]  has  generated  a  renewed  interest  in  this  direction. 

Another  relationship  has  been  proposed  by  Tokunaga  [45].  Using  Strahler’s  ordering 
system,  he  defined  e(_;  =  as  the  average  number  of  streams  of  order  j  flowing  laterally 
into  a  single  stream  of  order  i.  These  ratios  are  approximately  constant  independently  of 
order.  The  ratio  K=e(  /  is  also  approximately  constant. 

Using  elevations,  another  interesting  relationship  has  been  found  empirically  by  Flint 

[9],  The  relationship  between  slope  and  area  can  be  expressed  as  S  =  CA~e,  which  is,  like 

Horton’s  slope  law,  a  reflection  of  the  fact  that  stream  profiles  are  concave.  Gupta  and 
♦ 

Waymire  [11]  and  Tarboton  et  al.  [43]  have  also  studied  this  relationship.  Tarboton  et  al. 
[42]  has  postulated  a  relationship  of  the  form  Rb=R/2  based  on  fractal  arguments. 
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23.  Models  of  Channel  Network 

The  stability  of  the  Korton’s  ratios  under  a  variety  of  conditions  motivated 
geomorphologists  to  develop  random  models  that  reproduce  this  behavior.  Based  on  the 
work  on  topological  trees  by  Cayley  [7],  Shreve  [32]  developed  the  random  topology 
model.  In  this  modei  it  is  assumed  that  in  the  absence  of  geologic  controls,  channel 
networks  are  topologically  random,  i.e.  ail  the  topologically  distinct  channel  networks  with 
the  same  number  of  sources  are  equally  likely.  As  a  consequence  of  this  postulate,  in  such  a 
population  of  channel  networks,  the  most  probable  set  of  stream  numbers  Nu  approximately 
obeys  Horton’s  law. 

The  random  topology  model  was  extended  by  Smart  [36]  and  Shreve  [33],  [34]  to 
include  the  hypothesis  that  interior  and  exterior  link  lengths  as  well  as  their  associated  areas 
are  random  variables  with  separate  statistical  distributions.  This  new  postulate  makes  the 
model  have  very  interesting  properties  like  numerical  values  of  the  bifurcation,  length,  and 
area  ratios  that  are  similar  to  the  values  found  in  real  networks. 

Years  later.  Mesa  [29]  modelled  the  topologically  random  networks  as  a  birth  and 
death  Markov  process  where  at  any  distance  there  is  a  probability  of  bifurcation  or  death  of 
the  channel. 

All  the  properties  and  all  the  network  representations  which  have  been  shown  up  to 
this  point  have  a  common  characteristic:  they  look  at  a  snapshot  within  the  time  of  the 
network  history.  One  of  the  first  attempts  to  model  the  process  of  growth  of  a  network  was 
made  by  Howard  [15]  (see  also  van  der  Tak  [46]).  In  this  model  the  network  grows  over  a 
grid  of  points.  Nodes  are  defined  as  active  or  inactive.  Beginning  at  an  outlet  node  the 
network  extends  according  to  a  growth  probability  assigned  at  each  site.  Avoiding  closed 
loops,  the  algorithm  proceeds  until  it  achieves  a  certain  drainage  density  level. 

Another  model  is  the  random  walk  model  proposed  by  Leopold  and  Langbein 
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[21]  where  source  points  and  source  areas  assigned  to  them  are  distributed  as  a  Poisson 
process  in  space.  Then,  random  walks  proceed  from  the  sources  until  they  hit  the  boundary 
or  find  other  realizations  originating  in  other  sources.  In  this  way  networks  are  generated. 

Other  works  like  the  ones  by  Kirkby  [18]  and  Smith  and  Bretherton  [40]  look  at 
sediment  transport,  especially  at  the  hillslope  scale.  The  idea  is  to  consider  a  sediment 
continuity  equation,  with  flow  in  the  steepest  slope  direction.  The  amount  of  sediment 
transported  is  a  function  of  slope  and  flow  of  water.  The  amount  of  water  flowing  depends 
on  the  contributing  area  to  each  point.  Smith  and  Bretherton  [40]  present  a  perturbation 
analysis  of  the  2D  case  around  the  ID  equilibrium  solution.  They  find  that  small 
perturbations  grow  in  concave  2D  landscapes,  an  instability  that  is  manifested  in  channel 
growth.  Luke  [22]  generalizes  this  approach  to  include  erosion  and  deposition.  In  this  case 
the  instabilities  develop  into  shock  discontinuities  interpreted  as  channels. 

Recently  Willgoose  et  al.  [47]  have  developed  a  catchment  and  channel  network 
evolution  model.  This  model  is  used  in  this  work  and  will  be  explained  in  Section  2.4. 

2  J  A  Simulation  Model  of  Leaf  Veins  Growth 

Meinhardt  [24],  [25]  developed  a  model  that  simulates  the  growth  of  leaf  vein  cells 
and  their  differentiation  from  normal  leaf  cells.  Four  variables  were  used  in  the  model,  an 
activator  variable  "a”  that  triggers  the  differentiation  process;  an  inhibitor  h  that  has  a 
countereffect  and  reduces  production  of  activator;  a  substrate  z  consumed  by  the 
differentiation  process  and  a  differentiation  variable  Y  that  determines  whether  a  cell  is  a 
leaf  vein  cell  or  a  normal  cell.  The  model  studied  by  Meinhardt  was: 

dt  h  ax1 
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—  =  ca2z-  vh  +  D,j^  +  pxY  (^.6) 

Bt  Bx2 


d z  Kv  n  B2z 

--Co-ycvZ-SYz  +  D  — 


(2.7) 


—  =  c,a-OAY  + 
Bt  ' 


Y2 

1+9  Y2 


(2.8) 


In  the  simulation  process  the  initial  values  of  the  activator  and  inhibitor  fields  are 
small  and  z=l  everywhere.  A  point  is  chosen  as  the  seed  and  a  value  of  Y=1  is  defined  at 
this  point.  A  random  field  with  a  small  coefficient  of  variation  is  used  as  input  in  the 
parameter  c.  Willgoose  et  al.  [47]  studied  this  model  and  showed  analitically  how  very 
small  differences  in  the  c  field  generated  completely  different  leaf  vein  network  structures. 
The  explanation  was  found  in  the  chaotic  nature  of  the  system.  Chapter  4  will  explain  in 
more  detail  the  concept  of  chaos  in  a  dynamical  system.  Willgoose  et  al.  [47]  also  studied 
three  features  that  made  possible  the  generation  of  branched  structures  in  the  Meinhardt 
model:  a  potential  growth  region  exists  around  the  growing  branch  and  moves  with  it;  the 
growth  potential  is  supressed  behind  the  growing  tip  so  that  the  branch  grows  as  a  line;  a 
repulsion  between  growing  branches  exists  so  that  no  closed  loops  appear  and  the  final 
network  fills  the  entire  space. 


2.4  The  WBR  Model 

Based  on  the  conclusions  derived  in  the  study  of  the  Meinhardt  model,  Willgoose, 
Bras  and  Rodriguez-Iturbe  [47]  constructed  a  model  that  simulates  the  catchment  evolution 
and  network  growth  in  a  grid.  The  model  uses  some  concepts  of  erosion  engineering  to 
define  an  activator  function  that  triggers  the  growth  of  channels;  the  Einstein-Brown 
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sediment  transport  equation;  and  the  idea  of  preferred  erosion  in  channels.  Four  variables 
are  studied;  the  elevation  z,  the  indicator  Y  that  tells  whether  a  node  is  a  channel  or  a 
hillslope,  and  the  channel  formation  function,  a.  The  form  of  this  model,  henceforth  called 
the  WBR  model  is: 


^=vr,I,a/,+0_i 

(2.9) 

BY-  r  Y2  - 1 

(2.10) 

<*j  *  05  Qjml  Sjn 5 

(2.11) 

where  Q,  is  the  sediment  transport  (assumed  to  be  equal  to  (3t  Qjm\  Sf\  in  channels  and 
0fPi  Qj* t  Sf  1  with  0,<1  in  hillslopes);  Q;  is  the  overland  flow  which  is  proportional  to  the 
area  that  contributes  to  each  node;  S;  is  the  steepest  slope  downhill;  c0  the  tectonic  input  at 
node  j;  I-  an  indicator  function  that  says  whether  node  i  drains  into  node  j;  and  Kt,  Kj  and  d, 
are  constants.  The  physical  processes  on  which  this  model  is  based  will  be  explained  in 
more  detail  in  Section  3.3.  This  model  is  the  basis  for  the  subsurface  saturation 
modificadon  model  described  in  Chapter  3.  In  this  chapter  a  more  detailed  explanation  of 
the  different  variables,  the  physical  mechanisms  behind  them  and  details  about  the 
implementation  of  the  simulations  will  be  presented. 

The  WBR  model  generates  realistic  looking  catchments  and  channel  network  with 
values  of  statistical  measures  similar  to  the  networks  found  in  nature.  In  this  work  an 
extension  of  the  model  to  include  a  different  way  of  generating  overland  flow  and  a  closer 
look  at  the  chaotic  behavior  of  the  system  will  be  presented. 
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Chapter  3 

A  Catchment  Evolution  Model  with  a  Subsurface  Saturation  Mechanism 
3.1  Introduction 

In  this  chapter  a  modified  version  of  the  WBR  model  is  developed.  The  purpose  is  to 
modify  the  mechanism  of  overland  flow  production  to  include  subsurface  saturation  besides 
the  classical  Hortonian  infiltration  excess.  In  the  Hortonian  mechanism  overland  flow  is 
assumed  to  be  generated  uniformly  over  the  whole  area  whereas  in  the  subsurface 
saturation  mechanism  only  saturated  areas  contribute  to  the  overland  flow.  The  modified 
model  will  permit  us  to  understand  the  role  of  overland  flow  distribution  in  the  overall 
behavior  of  the  system. 

The  organization  of  this  chapter  is  as  follows.  An  analytical  formulation  of  the 
subsurface  saturation  mechanism  is  derived.  This  formulation  is  then  incorporated  in  the 
WBR  model  An  example  of  a  typical  simulation  and  the  variation  in  time  of  the  spatial 
distribution  of  the  principal  variables  of  the  model  are  shown.  Afterwards,  four  different 
random  elevation  fields  are  used  as  initial  conditions  for  different  smuiations  with  exactly 
the  same  parameter  values.  The  variability  of  common  geomorphological  statistics  is  then 

illustrated.  The  question  of  whether  chaos  is  the  explanation  for  this  variability  is  addressed 
in  Chapter  4. 

The  following  section  studies  the  influence  of  diffusive  sediment  transport  processes 
like  rainsplash  and  rockslide  in  the  evolution  of  the  networks.  Then  a  non-dimensional 
formulation  of  the  modified  model  is  developed.  In  this  formulation  the  original  parameters 
are  transformed  into  new  non-dimensional  ones  that  will  allow  comparisons  between 
catchments  with  different  scales. 
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whether  a  hillslope  node  is  saturated  or  not.  Then,  various  characteristics  of  importance  that 
have  been  studied  and  measured  in  real  basins  are  examined  in  the  resulting  networks  that 
are  generated  by  the  modified  version  of  the  WBR  model.  Finally,  in  Section  3.11  the 
influence  of  a  probabilistic  distribution  of  rainfall  on  the  overall  behavior  of  the  system  is 
studied. 

3.2  Overland  Flow  Distribution  Using  the  Subsurface  Saturation  Mechanism 

This  section  presents  the  analytical  formulation  of  the  procedure  to  assign 
contributing  flow  to  each  node  using  the  subsurface  saturation  mechanism.  The  overall 
reasoning  is  as  follows.  Assuming  an  exponential  profile  of  hydraulic  conductivity  with 
depth,  the  amount  of  water  required  to  obtain  saturation  at  each  node  can  be  calculated 
based  on  an  average  net  rainfall.  Then,  for  each  node  the  difference  between  the  rainfall  at 
the  mean  flood  event  and  the  saturation  deficit  can  be  calculated.  If  the  difference  is 
positive  then  overland  flow  is  produced  at  that  node.  By  adding  this  difference  along 
elevation  gradients  contributing  overland  flows  to  each  node  can  be  found. 

The  following  derivation  is  based  on  [4],  [5],  [6],  [8],  [19],  [20],  [30],  and  [35].  First 
of  all,  assume  exponential  decrease  of  hydraulic  conductivity  in  the  soil  profile  [2]: 

K(z)  =  K0e~a\z  (3.1) 

where  z  is  the  depth  from  the  soil  surface,  K(z)  and  K0  are  hydraulic  conductivities  at  depth 
z  and  at  the  soil  surface  respectively  and  ctj  is  a  parameter. 

If  the  water  table  depth  is  h  and  the  maximum  depth  of  interest  is  H,  then  the 
subsurface  flow  rate  at  node  i  is: 

qt  =  S^HK{z)dz  =  S,  [e~ai* - e“iw]  (3.2) 
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where  S;  is  the  steepest  downhill  slope.  Equation  (3.2)  assumes  that  the  local  hydraulic 
gradient  is  everywhere  equal  to  the  surface  slope  angle  [30].  If  H  is  large  enough,  then  q; 
can  be  approximated  to: 


(3.3) 


The  second  main  assumption  is  to  consider  a  moisture  profile  as  shown  in  Figure  3-1 
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Figure  3-1:  Soil  Profile  Assumption  in  Subsurface  Saturation  Mechanism 


The  saturation  deficit  is  defined  as: 

D^hie^-4 o) 


(3.4) 


Replacing  in  Equation  (3.3): 


r  An  r  i 

*  =  H*)exph^ 


(3.5) 


Define: 


„  *0  a  9--®0 

02  s;  p‘— sr 


(3.6) 


The  values  of  these  two  constants  are  assumed  to  be  homogeneous  over  the  whole 


basin.  Then  replacing  (3.6)  in  (3.5)  results  in: 
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(3.7) 


if  this  quantity  is  positive  or  D,-=0  otherwise.  A  negative  value  of  the  expression  (3.10) 
(i.e.D,-=0)  implies  that  the  node  i  is  saturated.  The  saturation  criteria  for  a  node  is  then: 


(3.11) 


For  every  point  the  excess  water  is  calculated  as  the  difference  between  P3  and  D,. 
is  considered  as  the  rainfall  at  the  mean  flood  event  (for  a  definition  of  this  concept  and  its 
relationship  with  the  sediment  transport  equation  to  be  used  in  the  model,  see  Appendix 
C.1.3.  in  [47]).  The  sum  along  the  elevation  gradient  of  the  excess  water  at  each  node  gives 
the  contributing  overland  flow  to  be  used  in  the  modified  modeL 
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3.3  Model  Description 


This  section  presents  the  subsurface  saturation  mechanism  modification  to  the  WBR 
model.  Two  variables  are  solved  on  a  plane  grid:  the  elevation  at  each  node  zt  and  the 
indicator  function  Y-  that  identifies  a  node  as  either  a  channel  or  a  hillslope.  The  downhill 
slopes  S i  determine  flow  directions.  Following  these  flow  directions  the  contributing  area 
A-  to  each  node  is  calculated.  The  next  step  is  to  find  the  deficits  D(  using  areas  and  slopes 
and  the  saturation  criteria  decribed  in  Section  3.2.  The  excess  water  is  summed  downhill 
along  the  flow  directions  giving  as  a  result  the  total  overland  flow  through  each  node  Q(. 
The  flows  Q,  enter  the  channel  formation  function  a,  and  sediment  transport  Q*,  equations. 
These  two  latter  variables  go  into  the  elevation  and  indicator  function  differential  equations. 


The  resulting  model  is: 


dfi 

dr 


+ 


1 

p,(i-")V 


(3.12) 
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(3.14) 

(3.15) 

(3.16) 


(3.17) 


-24- 


T 


(3.18) 


where  <x>  stands  for  max(x,0)  and  the  subindices  i  and  j  denote  grid  nodes.  Most  of  the 
variables  have  already  been  defined.  The  following  is  a  complete  list  of  all  constants  and 
variables  involved  in  the  equations: 


z(  =elevation 

Y(  =indicator  variable  for  channelization  (Ohillslope,  1  channel) 

Qs  =sediment  transport 

a,  =channel  formation  function  function 

Q,  =discharge 

D(  =saturation  deficit 

t  =time 

C0  =tectonic  input 

i 

=density  of  eroded  material 

n  =porosity  of  material  before  erosion  and  after  deposition 
L?  =grid  spacing 

=indicator  function  for  whether  node  j  drains  into  or  out  node  i 
-l(j=i,drainage  of  node  i),  l(node  j  drains  into  node  i), 

0(node  j  does  not  drain  into  or  out  of  node  i) 

D:  =diffusivity  constant  for  diffusive  transport  processes 

\j  horizontal  distance  in  direction  j 

d,  =rate  constant  for  channel  growth 

c,  =l/threshold  on  the  channel  formation  function  of  channelization 

m,,n,  =powers  of  Q-  and  S(  in  the  sediment  transpon  equation 
P5  =multiplicative  factor  on  channel  formation  function  equation 

m5,n5  =powers  of  Q(  and  S,  in  the  channel  formation  function  equation 
P3  =additive  term  on  overland  flow  equation 

P6  =multiplicative  factor  on  saturation  deficit  equation 

Tc  =threshold  value  in  saturation  criteria 

Let  us  explain  in  more  detail  the  equations  that  govern  the  evolution  of  the  system. 
There  are  three  basic  variables:  elevation,  channelization  ind  cator  and  value  of  channel 
formation  function.  The  first  two  are  the  variables  to  be  solved  in  the  grid  plane.  They  are 
related  to  each  other  through  the  channel  formation  function  function.  Grid  points  are 
defined  as  channel  or  hillslope  nodes  according  to  the  value  of  the  indicator  function  Y(  (1 
for  channels,  0  for  hillslopes).  The  differential  equation  (3.13)  for  Y(  is  constructed  in  such 
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a  way  that  the  values  0  and  1  are  attractive  fixed  points.  Y(  depends  on  the  values  of  the 
channel  formation  function  function  defined  by  Equation  (3.15).  The  purpose  of  the  channel 
formation  function  is  to  trigger  the  process  that  transforms  a  hillslope  node  into  a  channel 
node.  Once  the  channel  formation  function  exceeds  a  threshold  value  given  by  1/c,,  the 
indicator  function  moves  from  0  to  1  at  a  rate  determined  by  d(.  The  specific  way  in  which 
Equation  (3.13)  rules  this  transformation  is  not  important  for  the  overall  behavior  of  the 
system:  any  equation  with  0  and  1  as  attractive  fixed  points  and  a  dependence  on  a 
parameter  that  behaves  like  a  will  work  as  well. 

The  definition  of  the  channel  formation  function,  in  Equation  (3.15),  is  a 
conceptualization  of  erosion  processes.  Concepts  like  overland  flow  velocity  and  threshold 
bed  shear  stress,  that  are  commonly  used  in  sedimentation  engineering,  can  be  placed  into 
the  framework  of  Equation  (3. 15)  (see  [47]  for  details). 

On  the  other  hand,  the  variation  in  grid  node  elevations  is  produced  by  three  factors: 
the  tectonic  input,  the  fluvial  sediment  transport  and  diffusive  transport.  The  tectonic  input 
can  take  a  variety  of  forms.  The  expressions  studied  in  this  work  correspond  to  a  single 
uplift  event  (see  [26])  or  a  continuous  process. 

The  fluvial  sediment  transport  term  is  simply  a  sediment  continuity  equation.  The 
amount  of  sediment  that  is  washed  away  from  the  node  by  fluvial  transport  processes  is 
calculated  using  Equation  (3.14).  The  term  f(Yt)  in  this  equation  parametrizes  the  difference 
in  sediment  transport  phenomena  that  occurs  in  hillslopes  and  channels.  A  small  value  of  O, 
in  Equation  (3.17)  describes  the  fact  that  sediment  transport  in  hillslopes  is  smaller  than  in 
channels.  Details  about  the  derivation  of  Equation  (3.14)  from  the  well  known  Einstein- 
Brown  sediment  transport  equation  can  be  found  in  [47], 

The  term  Q(  in  Equation  (3.14)  is  the  net  overland  flow  contributing  to  node  i. 
Equation  (3.16)  is  the  overland  flow  continuity  equation,  where  flow  from  adjacent  nodes 
that  drain  into  node  i  and  excess  water  from  that  node  are  equated  to  the  flow  from  the 
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node.  The  excess  water  at  each  node  is  calculated  based  on  the  saturation  deficit  equation 

(3.18)  as  described  in  Section  3.2. 

The  last  term  in  the  differential  equation  (3.12)  for  the  elevation  field  corresponds  to 
transport  processes  that  occur  on  the  hiilslope  and  can  be  modelled  as  diffusive  processes. 
These  processes  are  soil  creep,  rainsplash  and  rockslide. 

Finally,  the  sediment  transport  is  the  non-linear  link  between  the  elevation  and 
indicator  variables.  The  sediment  movement  changes  node  elevations  and  this  affects  both 
slopes  and  contributing  overland  flow.  Then  the  channel  formation  function  changes  and  its 
new  value  enters  the  indicator  differential  equation,  triggering  the  channelization  process  if 
the  value  surpasses  the  threshold  1/c,.  The  indicator  function  is  not  just  an  artifact  to 
differentiate  hiilslope  and  channel  nodes  but  it  also  affects  the  way  sediment  transport  takes 
place.  The  preferred  transport  in  channels  creates  valleys  around  the  channels  and  the  whole 
slope  field  behaves  accordingly.  The  dependence  cycle  in  the  system  is  then  closed. 

In  the  implementation  of  the  program,  the  system  is  solved  in  a  rectangular  grid.  An 
initial  random  elevation  field  with  a  very  small  coefficient  of  variation  is  generated  and  the 
elevation  at  the  bottom  left-hand  comer  is  set  at  a  lower  value.  That  node  is  defined  as  the 
outlet  and  its  elevation  is  kept  constant  during  the  simulation. 

3.4  A  Simulation  Example 

In  order  to  show  the  overall  behavior  of  the  system  described  by  Equations  (3.12)  to 

(3.18) ,  a  simulation  example  will  be  presented  in  this  section.  No  detailed  analysis  of  the 
different  network  properties  will  be  performed  here.  Such  an  analysis  will  be  presented  in 
subsequent  sections.  Rather,  the  variation  in  time  of  the  most  important  properties  of  the 
catchment  and  the  channel  network  for  a  typical  s  .ulation  are  going  to  be  shown 
graphically. 
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The  parameters  for  this  simulation  BR643SATU10  are  shown  in  Appendix  A.  Four 
steps  in  the  evolution  of  the  network  have  been  chosen  to  show  the  variation  in  time  of 
different  properties.  The  time  interval  spans  over  the  whole  network  growth  period.  The 
properties  to  be  examined  are  elevation,  slope,  contributing  area,  overland  flow  and  channel 
formation  function  at  each  node.  For  a  better  appreciation  of  the  properties,  three  things  are 
presented  in  each  graph:  the  spatial  distribution  of  the  property  underneath  an  isometric 
view  of  the  elevation  field:  a  contour  map  corresponding  to  the  values  of  the  property  under 
study;  and  the  simulated  channel  network  at  that  time  step.  The  contour  levels  in  each 
figure  are  not  evenly  spaced:  instead,  their  values  have  been  chosen  to  highlight  certain 
aspects  of  the  behavior  of  the  system. 

Figures  3-2  to  3-5  show  the  evolution  of  the  elevation  field.  The  preferential  erosion 
at  channels,  determined  by  O,  ,  produces  well  differentiated  valleys  along  the  channel. 
Figure  3-6  to  3-9  correspond  to  the  steepest  slope  S,  in  the  downhill  direction  at  each  node. 
Notice  that  the  steepest  slopes  in  the  field  occur  not  on  the  channel  heads  but  on  the  lateral 
hills  along  the  channel.  The  zone  outside  the  region  captured  by  the  network  has  slopes  that 
are  very  small.  Figures  3-10  to  3-13  present  the  contributing  area  Ar  As  expected,  channel 
nodes  have  the  larger  areas  as  a  result  of  the  organization  of  flow  directions  in  the 
catchment.  Figures  3-14  to  3-17  show  which  nodes  are  saturated  (in  black)  and  which 
nodes  are  unsaturated  (in  white).  Equation  (3.1 1)  was  the  criteria  for  selecting  the  saturation 
status  of  each  node.  At  the  beginning,  every  node  was  saturated:  with  time,  only  channel 
nodes  and  nodes  around  them  remain  saturated.  Figures  3-18  to  3-21  correspond  to  the 
channel  formation  function  distribution.  The  contours  show  how  the  higher  values  of  the 
channel  formation  function  are  located  in  those  zones  of  active  network  growing  at  the 
channel  heads.  In  the  original  WBR  model  the  areas  and  contributing  flow  were  directly 
proportional.  In  the  subsurface  saturation  mechanism  modification  the  overland  flow  is  not 
directly  proportional  to  areas  anymore,  even  though  strong  links  exist  between  both 
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variables.  In  the  modified  model  overland  flows  are  calculated  according  to  the  derivation 
presented  in  Section  3.2.  Finally,  Figure  3-22  shows  the  variation  in  time  of  the 
hypsometric  curve  for  the  catchment.  Detailed  analysis  of  this  curve  under  different 
parameter  values  will  be  presented  in  a  later  section. 

Table  3-1  presents  the  values  of  commonly  studied  statistics  for  the  final  network. 
Notice  that  all  the  values  of  the  statistics  are  within  reasonable  limits. 

The  overall  behavior  of  the  system  variables  have  been  shown  in  a  typical  simulation. 
The  final  result  is  a  realistic  catchment  and  geomorphologic  network  with  characteristics 
similar  to  real  networks.  An  analysis  of  the  dependence  of  the  system’s  behavior  on  the 
parameters  and  a  comparison  between  the  original  WBR  model  with  its  subsurface 
satuation  modification  will  be  presented  in  the  following  sections. 

3-5  Variability  in  Final  Channel  Networks 

This  section  will  illustrate  the  high  variability  present  in  the  final  results  of  the 
catchment  evolution  model.  This  variability  is  the  reason  for  an  apparent  randomness  in  the 
final  channel  network.  The  differential  equations  of  the  model  have  two  important 
characteristics:  they  are  non-linear  and  they  are  coupled  in  space.  This  is  the  reason  why 
even  if  two  different  simulations  have  the  same  parameter  values,  small  differences  in  the 
elevation  field  generate  completely  different  networks. 

Four  simulations  are  shown  in  this  section.  All  of  them  have  the  same  parameter 
values.  The  difference  between  simulations  lies  in  the  initial  elevation  field.  Different 
realizations  of  random  elevation  fields  with  identical  statistical  properties  were  used.  Final 
networks  are  shown  in  Figure  3-23. 

For  each  of  the  four  catchment  evolutions  common  geomorphological  measures  were 
examined:  Strahier’s  bifurcation,  length,  slope  and  area  ratios.  The  evolution  in  time  of 


Figure  3-4:  Elevation  Spatial  Distribution  Timestep  5000 
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Figure  3-5:  Elevation  Spatial  Distribution  Timestep  7000 


Figure  3. 


Figure  3-7 :  Slopes  Spatial  Distribution  Timestep  3000 


Figure  3-8:  Slopes  Spatial  Distribution  Timestep  5000 
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Figure  3-9:  Slopes  Spatial  Distribution  Timestep  7000 


Figure  3-12:  Contributing  Area  Spatial  Distribution  Timestep  5000 
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Figure  3-13:  Contributing  Area  Spatial  Distribution  Timestep  7000 
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Figure  3-16:  Saturated  and  Unsaturated  Nodes  Timestep  5000 
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Figure  3-17:  Saturated  and  Unsaturated  Nodes  Timestep  7000 
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Figure  3-18:  Channel  Formation  Function  Spatial  Distribution  Timestep  1000 
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Figure  3-19:  Channel  Formation  Function  Spatial  Distribution  Timestep  3000 
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Figure  3-20:  Channel  Formation  Function  Spatial  Distribution  Timestep  5000 
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Figure  3-21:  Channel  Formation  Function  Spatial  Distribution  Timestep  7000 
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Figure  3-22:  Hypsometric  Curve  Variation  in  Time 


Statistic 

Value 

R*(l-2) 

5.08 

Rj(l-2) 

2.85 

Rl(1-2) 

2.60 

R,(l-2) 

5.85 

Rj,(2-3) 

4.00 

Rj(2-3) 

2.15 

Rt(2-3) 

1.82 

R/t(2-3) 

4.96 

R* 

4.75 

R, 

2.67 

Rl 

2.35 

R„ 

5.53 

K 

1.96 

2.67 

D/ 

11.54 

Magnitude 

61 

Mean  Relief 

9.90 

Mean  Stream  Relief 

3.12 

Table  3-1:  Final  Network  Statistic  Values  Simulation  BR643SATU10 
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Figure  3-23:  Final  Catchments  for  Simulations  with 
Different  Initial  Elevation  Fields 
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each  of  these  measures  appears  in  Figures  3-24  to  3-31.  These  ratios  were  calculated 
between  links  of  order  1  and  2  and  between  links  of  order  2  and  3  according  to  Strahler’s 
ordering  system.  Notice  that  even  though  the  values  are  within  reasonable  limits,  a  scatter 
of  results  is  apparent. 

Figures  3-32,  3-33  and  3-34  correspond  to  the  time  evolution  of  drainage  density, 
magnitude  and  mean  link  lengths  respectively.  Notice  how  the  variability  in  drainage 
density  is  much  smaller  than  the  variability  of  those  measures  based  on  topological 
variables. 

Finally,  Figure  3-35  shows  the  variation  of  the  hypsometric  curve  with  time.  The 
variation  of  both  the  hypsometric  curves  and  the  drainage  density  for  the  different 
evolutions  is  relatively  small.  The  advantage  of  these  two  measures  is  that  they  are 
physically  rather  than  topologically  based. 

Chapter  4  will  look  in  greater  detail  at  this  inherent  variability  of  the  catchment 
evolution  model.  This  variability  will  be  related  to  the  concept  of  transient  chaos 
quantitative  measures  will  illustrate  this  effect. 

3.6  Hypsometric  Curves  for  Cases  with  Continuous  and  Instantaneous  Uplift 

In  this  section  the  time  evolution  of  the  hypsometric  curves  for  simulations 
BR643SATU10  and  BR643SATU33  are  studied.  Both  simulations  have  the  same 
parameter  values.  The  difference  between  them  is  in  the  way  tectonic  uplift  is  applied.  In 
simulation  BR643SATU10  the  uplift  is  given  as  a  single  event  at  the  beginning  of  the 
evolution  while  in  simulation  BR643SATU33  the  tectonic  input  is  a  continuous,  constant 
and  uniformly  distributed  process  in  time.  Figures  3-36  and  3-37  present  the  hypsometric 
curve  evolutions  for  both  cases.  The  jump  at  the  upper  part  of  the  curves  is  an  effect  of  the 
lower  outlet  elevation.  With  time,  the  elevation  jump  at  the  outlet  comer  is  reduced  because 
of  the  sediment  transport,  and  a  smoother  field  of  elevations  appears. 


Ratio.  Order  1-2 


Figure  3-25:  Strahler  Bifurcation  Ratio.  Order  2-3 
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Figure  3-27:  Strahier  Length  Ratio.  Order  2-3 
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Figure  3-29:  Strahler  Slope  Ratio.  Order  2-3 


Figure  3-31:  Strahier  Area  Ratio.  Order  2-3 


Figure  3*32:  Time  Evolution  of  Drainage  Density 


Figure  3-33:  Time  Evolution  of  Magnitude 


1  2  3  4  .  5  6  7  8 

Time  103 


Figure  3-34:  Time  Evolution  of  Mean  Link  Length 
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Figure  3-35:  Time  Evolution  of  Hypsometric  Curve 
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Figure  3-36:  Hypsometric  Curve  Evolution. 
Episodic  Uplift  Case  10 


Figure  3-37:  Hypsometric  Curve  Evolution. 
Continuous  Uplift  Case  33 
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The  case  with  continuous  tectonic  uplift  tends  to  evolve  into  an  S-shaped 
hypsometric  curve.  This  form  appears  because  the  tectonic  input  is  applied  everywhere  in 
the  elevation  field  at  the  same  rate.  The  erosion  at  hillslopes  with  high  elevation  is  small 
because  the  subsurface  saturation  threshold  reduces  contributing  overland  flow.  Therefore, 
those  nodes  with  higher  elevation  in  the  catchment  increase  in  their  elevation  while  lower 
nodes  decrease  further  in  elevation  to  the  point  of  dynamic  equilibrium  where  the  S-shape 
is  attained.  Notice,  however,  that  the  variable  shown  is  normalized  elevation. 

The  case  with  instantaneous  tectonic  uplift  has  a  kink  that  is  created  because  of  the 
differences  in  sediment  transport  between  saturated  and  unsaturated  nodes.  It  will  be  seen  in 
the  next  section  how  diffusive  sediment  transport  processes  have  a  major  influence  on  the 
evolution  of  the  hypsometric  curves. 

3.7  Influence  of  Diffusion  in  Simulations 

The  evolution  equation  (3.12)  for  elevation  in  the  catchment  evolution  model  has 
three  components  of  which  the  two  most  important  are  the  tectonic  input  and  the  fluvial 
sediment  transport.  The  third  component  corresponds  to  diffusive  processes  like  rainsplash 
and  rockslide.  Their  effect  is  specially  important  on  unsaturated  hillslope  nodes.  While  in 
the  original  WBR  model  overland  flow  occurred  at  every  node,  in  its  subsurface  saturation 
modification  there  is  no  overland  flow  in  unsaturated  nodes.  It  is  at  those  nodes  where 
diffusive  processes  are  the  impotant  mechanism  for  sediment  transport. 

Simulations  BR643SATU10  and  BR643SATU32  share  the  same  parameter  values 
except  for  the  value  of  the  diffusion  coefficient  Dr  In  the  first  case  Dz=0  while  in  the 
second  case  Dr=3.5xlO-5.  Figure  3-38  shows  the  temporal  evolution  of  the  drainage  density 
for  both  simulations. 


Figure  3-39  corresponds  to  the  hypsometric  curves  as  time  varies.  Notice  how  in  the 
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case  with  diffusion,  hypsometric  curves  decrease  faster  because  of  the  increased  sediment 
transport  that  comes  especially  from  unsaturated  nodes.  In  contrast  with  the  original  WBR 
model  where  overland  flow  was  present  at  every  node  and  therefore  the  influence  of 
diffusion  was  small,  in  the  subsurface  saturation  modification  the  diffusive  sediment 
transport  processes  have  an  interesting  influence  on  the  evolution  of  hypsometric  curves.  In 
those  cases  with  Dz=0  nodes  at  high  elevation  do  not  erode  as  fast  -or  do  not  erode  at  all-  as 
lower  nodes  that  usually  have  larger  contributing  areas  and  smaller  slopes  which  makes 
them  susceptible  to  saturation.  The  presence  of  diffusion  makes  the  difference  in  behavior 
between  both  sets  of  nodes  smaller,  smoothing  out  the  kink  that  appears  in  the  cases  with 
no  diffusion. 


Despite  the  influence  that  diffusive  processes  have  on  the  evolution  of  hypsometric 
curves,  physically-based  geomorphological  measures  do  not  differ  very  much  in  both  cases, 
which  reinforces  the  preference  of  these  measures  over  the  topologically  based.  Notice, 
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however,  that  the  form  of  the  final  networks  is  completely  different;  it  also  results  in 
different  values  of  topologically-based  network  measures.  Figure  3-40  presents  the  form  of 
the  final  networks.  Chapter  4  will  examine  how  small  perturbations,  in  this  case  produced 
by  a  new  sediment  transport  factor,  can  affect  the  entire  catchment  evolution. 

3.8  Non-Dimensional  Formulation  of  the  Model 

Given  that  the  different  parameters  and  variables  of  the  model  are  strongly 
interconnected  in  the  differential  equations,  in  order  to  make  thoughtful  comparisons 
between  basin  evolutions  in  different  points  of  the  parameter  space,  a  non-dimensional 
version  of  the  model  has  to  be  defined.  The  basic  idea  is  to  combine  parameters  that  work 
in  the  same  phenomena  and  to  generate  lumped  non-dimensional  parameters.  This  section 
is  mostly  based  on  [47]. 
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Figure  3-40:  Final  Networks. 

Cases  10  and  32 

The  first  step  is  to  define  the  basic  non-dimensional  variables: 


(3.19) 

(3.20) 

(321) 


(3.22) 
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where  z',x',R' and  t’  are  non-dimensional  elevation,  horizontal  distance,  runoff  rate  and  time 
respectively.  L„  Lx,  L^,  TR  and  T  are  the  corresponding  scales  in  elevation,  horizontal 
length,  runoff  depth,  runoff  time  and  catchment  evolution  time  respectively.  The 
expressions  are  self  explanatory,  except  perhaps  Equation  (3.21)  that  comes  from  the 
assumption  of  uniform  rainfall  over  the  area,  which  in  terms  of  the  model  implies  R=P3. 


The  second  step  is  to  define  second-level  variables  of  the  model  in  this  framework: 
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where  S',A',Qs',a,f'(Y)  and  Tc'  are  the  nondimensional  expressions  of  the  analogous 
variables  in  the  model.  The  only  factors  in  the  above  expressions  that  remain  to  be  defined 
are: 


P/(-t,y)  = 


P,UoO 


(3.30) 
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ct2(x,y) 
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(3.32) 


C0'(x,y)  = 


C0(x,y) 


(3.33) 


The  Equations  (3.30)  to  (3.33)  express  the  fact  that  P,,  p5,  a:  and  C0  may  vary  in 
space.  Therefore,  they  are  defined  in  terms  of  a  mean  scale  for  each  variable.  This  scale 
cannot  be  related  to  the  basic  non-dimensional  variables  because  they  correspond  to 
different  processes:  P,  and  P5  to  details  of  the  sediment  transport  phenomena  like  particle 
size;  c0  to  the  tectonic  uplift  process;  and  a,  to  hydraulic  conductivity  decrease  in  the  soil. 

The  last  step  in  the  process  is  to  lump  the  original  parameters  in  order  to  form  the 
non-dimensional  parameters.  The  way  the  original  parameters  are  combined  comes  from 
the  specific  form  of  the  equations.  The  non-dimensionalized  governing  equations  are  now: 
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(3.35) 


which  together  with  the  non-dimensional  subsurface  saturation  criteria  and  the  definitions 
of  second-level  variables  (QaQ^.  and  D)  forms  the  model.  The  saturation  criteria  in  non- 
dimensional  form  is: 


TSC  7/  <  j 


(3.36) 


TT,TS,  TD,  TC,  TA  and  TSC  are  the  nondimensional  tectonic  uplift,  sediment 
transport,  diffusion,  channelization  rate,  channel  formation  function  and  saturation  criteria 
numbers  respectively.  Their  definitions  follow  from  the  model  equations  and  are  now 
presented: 
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TSC  = 


-rLz 


(3.42) 


Based  on  these  non-dimensional  numbers,  comparisons  can  be  made  between 
catchment  evolutions.  When  two  catchments  have  identical  nondimensional  parameters  as  t 

well  as  identical  distributions  of  nondimensional  catchment  properties  like  discharge  Q\ 
sediment  transport  f'(Y),  elevation  z  and  tectonics  C'„,  they  can  be  considered  physically 
similar. 

Willgoose  et  al.  [ 47 1  present  two  definitions  of  non-dimensional  numbers.  The 
definitions  presented  before  correspond  to  transient  conditions.  There  are  analogous 
expressions  for  catchments  in  dynamic  equilibrium. 


3.9  Influence  of  TSC  in  the  Catchment  Evolution 

The  value  of  the  threshold  that  determines  whether  a  node  is  saturated  or  not  has 
important  implications  on  the  overall  behavior  of  the  system.  The  value  of  TC  (and  in  a 
broader  sense  the  value  of  TSC)  affects  the  amount  of  overland  flow  that  is  distributed  to 
the  nodes.  In  this  section  three  different  simulations  with  varying  values  of  TSC  will  be 
shown.  The  cases  to  be  presented  are  BR643SATU10,  38  and  39  which  have  values  of  TSC 
equal  to  0.92\103,  4.6xl03  and  0.28xl03  respectively.  All  the  other  parameter  values  are  the 
same.  Figures  3-41  to  3-43  show  an  isometric  view  of  the  catchment  once  the  network  stops 
growing,  a  contour  of  elevations  and  the  final  network. 

With  decreasing  values  of  TSC,  a  larger  number  of  grid  nodes  enter  into  the  saturated 
group  where  overland  flow  and  fluvial  sediment  transport  is  produced.  An  increase  of 
overland  flow  production  produces  a  higher  value  of  channel  formation  function,  triggering 
the  growth  of  new  channels  as  can  be  seen  by  comparing  Figures  3-42  and  3-41.  Channel 
valleys  are  wider  in  case  39  compared  with  case  10  because  of  the  increased  fluvial 
transport  in  hillslope  nodes  surrounding  channels. 
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If  the  value  of  TSC  is  small  enough,  then  all  the  nodes  would  be  saturated  and  the 
behavior  of  the  modified  model  would  be  identical  to  the  original  WBR  model.  Finally,  it 
should  be  noticed  that  it  is  the  difference  in  overland  flow  distribution,  produced  by  the 
saturation  deficits,  that  makes  analytical  calculations  more  difficult  in  the  subsurface 
saturation  case  compared  to  the  original  model. 

3.10  Analysis  of  Geomorphoiogical  Relationships  in  Simulated  Catchments 

In  this  section  a  series  of  hypotheses  and  relationships  that  have  been  studied  and 
tested  in  real  basins  will  be  examined  in  the  final  catchments  generated  by  the  subsurface 
modification  to  the  WBR  model. 

The  relationships  studied  are:  the  hypothesis  of  area-slope  renormalization  in 
channels  and  the  possibility  of  using  analysis  of  this  kind  to  distinguish  channels  from 
hillslopes  in  digital  elevation  maps;  statistical  properties  of  link  lengths  and  contributing 
areas;  and  the  behavior  of  some  fundamental  scales  of  the  catchment. 

3.10.1  Area-Slope  Renormalization 

Flint  [9]  found  empirically  that  slopes  and  channels  scale.  He  found  a  relationship 
between  area  A  and  slope  S  of  the  form: 

S  =  CA~®  (3.43) 

with  0  ranging  from  0.37  to  0.83  with  a  mean  of  0.6.  Recently,  Tarboton  [44]  also  found 
that  this  scaling  relationship  holds  between  areas  and  slopes  in  catchments  using  digital 
elevation  maps.  The  scaling  relation  was  tested  for  simulations  with  the  WBR  model  in 
[47]  and  an  excellent  fit  was  found  in  catchments  at  dynamic  equilibrium. 

In  this  subsection  simulations  using  the  subsurface  saturation  modification  are 
examined.  Figure  3-44  shows  the  data  points  of  area  and  slope  for  channelized  nodes.  The 
simulation  case  under  consideration  was  BR643SATU33  at  dynamic  equilibrium. 


Figure  3-41:  Elevation  Spatial  Distribution.  Case  10 
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<or 


Figure  3-42:  Elevation  Spatial  Distribution.  Case  38 


Figure  3-43:  Elevation  Spatial  Distribution.  Case  39 
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Area 


Figure  3-44:  Area-Slope  Renormalization. 

Case  33  at  Dynamic  Equilibrium 

The  fit  of  Equation  (3.43)  is  also  excellent  in  this  modified  model.  The  observed 
value  of  9  is  0.38.  The  analysis  of  Willgoose  et  al.  [47]  to  calculate  the  value  of  9  from  the 
parameters  of  the  model  cannot  be  applied  to  the  subsurface  saturation  modification.  The 
reason  is  that  in  the  modified  model  there  is  no  a  simple  and  direct  relationship  between  Q 
and  A  as  in  the  original  model  (Q=P3A)  because  of  the  presence  of  saturation  deficits. 

Tarbcton  [44]  hypothesizes  that  channels  and  hillslopes  can  be  distinguished  in  a 
digital  elevation  map  based  on  a  slope  break  that  appears  in  the  renormalization  analysis 
when  every  pixel  in  the  map  is  studied.  The  indicator  function  which  defines  exactly 
whether  a  node  is  channelized  or  not  in  the  WBR  model  permits  this  hypothesis  to  be 
tested.  Willgoose  et  al.  [47]  found  that  the  criteria  for  hillslope  distinction  was  not  a  break 
in  slope;  instead  he  found  a  separate  hillslope  renormalization  line.  This  behavior  of  the 
model  has  to  do  with  the  fact  that  the  only  difference  in  sediment  transport  mechanisms 
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between  channels  and  hillslopes  is  the  multiplicative  factor  Or  Figures  3-45  and  3-46  show 
the  area-slope  renormalization  analysis  tor  the  case  with  continuous  tecronic  uplift  before 
the  system  reaches  dynamic  equilibrium.  All  the  grid  nodes  are  included  in  these  graphs.  It 
can  be  seen  at  Figure  3-45  that  at  time  4000  transient  points  are  present  in  the  bottom  left- 
hand  comer.  These  points  correspond  to  those  nodes  where  the  action  of  the  network  has 
not  arrived.  As  time  passes,  the  dispersion  of  the  points  in  the  graph  reduces  until  it  reaches 
the  straight  line  renormalization  at  dynamic  equilibrium. 


Figure  3-45:  Area-Slope  Renormalization.  All  Nodes  Included. 

Case  33  Before  Dynamic  Equilibrium.Time  4000 

Figures  3-47,  3-48  and  3-49  show  the  area-slope  graphs  at  three  different  times  in  the 
evolution  ot  a  system  that  includes  diffusive  sediment  transport  processes  but  instantaneous 
tectonic  uplift.  Again,  there  are  transient  points  with  small  area  and  very  small  slope.  These 
points  in  the  graph  dissapear  later  in  the  evolution  from  the  effect  of  sediment  wash. 
However,  even  though  the  fit  of  expression  (3.43)  is  good  as  can  be  seen  in  Figure  3-49.  it 
is  not  as  precise  as  it  was  in  the  continuous  uplift  case.  In  the  diffusion  case,  a  reduction  of 
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slope  is  present  as  areas  decrease.  Later  on,  some  hillslope  nodes  are  left  behind  in  the 
sediment  wash  process  because  they  belong  to  the  unsaturated  group. 

3.10.2  Link  Lengths  of  Simulated  Catchments 

This  section  will  study  the  planar  properties  of  the  simulated  catchments  using  the 
subsurface  saturation  modification.  The  results  will  be  compared  to  the  simulations 
performed  with  the  original  model.  Also  some  relationships  that  have  been  studied  by  field 
workers  will  be  examined.  The  relationships  to  be  considered  are:  that  the  probability 
distribution  of  exterior  and  interior  link  lengths  are  different  [37];  that  the  link  length 
distribution  is  exponential  [37]  or  Gamma  [46];  and  that  link  lengths  are  correlated  with 
link  magnitude  or  the  downstream  link  magnitude  [39]  and  [1]. 

In  order  to  test  the  different  hypotheses,  Table  3-0  presents  the  mean  link  length  for 
exterior,  interior  and  the  whole  set  of  links  (denoted  by  Le,  L,  and  L,  respectively). 
Following  [47]  a  t-test  was  performed  in  order  to  find  out  whether  a  statistically  significant 
difference  existed  between  interior  and  exterior  mean  lengths.  The  null  hypothesis  is: 

LJLl  -  1 

H0  ;  - -  =  0  (3.44) 

®  L  /L 

e'i 

The  mean  value  of  the  test  variable  is  0.29.  Therefore  the  hypothesis  cannot  be 
rejected  at  the  5%  confidence  limit,  i.e.  no  statistical  difference  was  found  between  the 
mean  values  of  exterior  and  interior  link  lengths. 

Based  on  this  result,  it  is  possible  to  study  the  distribution  of  all  the  links  at  once, 
both  interior  and  exterior.  Both  the  mean  and  the  standard  deviation  of  L,  for  different 
simulations  are  presented  in  Table  3-III.  Also  in  this  table  the  shape  factor  for  a  moment- 
fitted  gamma  distribution  was  calculated.  The  mean  value  of  the  shape  factor  was  1.43  ± 
0.28  which  is  consistent  with  the  hypothesis  that  link  lengths  are  gamma  distributed  with  a 
shape  parameter  factor  of  about  1.5  to  2.0  [46]. 
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SIMULATION 

4 

L,. 

l_At l 

°LJLi 

BR643SATU10 

3.81 

3.14 

+0.85 

BR643SATU21 

2.71 

2.62 

+0.14 

BR643SATU29 

3.91 

3.83 

+0.08 

BR643SATU30 

2.63 

3.36 

-0.86 

BR643SATU31 

3.03 

2.49 

+0.86 

BR643SATU32 

3.42 

3.19 

+0.29 

BR643SATU33 

2.74 

3.10 

-0.46 

BR643SATU34 

4.59 

3.44 

+  1.33 

BR643SATU35 

4.20 

4.02 

+0.18 

BR643SATU36 

3.22 

3.37 

-0.18 

BR643SATU37 

4.36 

3.90 

+0.47 

BR643SATU39 

3.28 

3.22 

+0.07 

BR643SATU40 

5.09 

7.50 

-1.28 

BR643SATU42 

5.17 

2.81 

+3.34 

BR643SATU43 

2.58 

2.55 

+0.05 

BR643SATU44 

6.55 

5.52 

+0.74 

BR643SATU45 

2.21 

2.69 

-0.70 

Table  3-II:  T-test  for  Mean  Exterior  and  Interior  Link  Lengths 


Finally,  link  lengths  were  plotted  against  magnitude  of  the  same  link  and  the 
following  link  downstream.  Figures  3-50  and  3-51  correspond  to  these  plots.  The  case 
under  consideration  corresponds  to  runs  BR643SATU10,34,35  and  36  which  have  the  same 
different  parameters  but  different  initial  random  elevation  fields.  Correlations  for  case  10 
were  found  to  be  low  (-0.15  for  link  length  vs.  magnitude  and  -0.18  for  link  length  vs. 
downstream  magnitude)  in  agreement  with  results  by  [47]  and  [1]. 


SIMULATION 

L, 

BR643SATU10 

3.48 

3.04 

1.31 

BR643SATU21 

2.66 

2.17 

1.50 

BR643SATU29 

3.87 

3.21 

1.45 

BR643SATU30 

2.99 

2.33 

1.65 

BR643SATU31 

2.77 

2.20 

1.59 

BR643SATU32 

3.31 

3.21 

1.06 

BR643SATU33 

2.91 

2.11 

1.90 

BR643SATU34 

4.02 

3.39 

1.41 

BR643SATU35 

4.11 

3.88 

1.12 

BR643SATU36 

3.29 

2.34 

1.41 

BR643SATU37 

4.13 

3.44 

1.44 

BR643SATU39 

3.25 

2.68 

1.21 

BR643SATU40 

6.26 

5.89 

1.13 

BR643SATU42 

4.00 

4.27 

0.88 

BR643SATU43 

2.57 

1.98 

1.68 

BR643SATU44 

6.05 

4.36 

1.93 

BR643SATU45 

2.45 

1.89 

1.67 

Table  3-IH:  Shape  Parameter  of  a  Gamma  Distribution  fitted  to  Link  Lengths 

3.10.3  Link  Contributing  Areas  in  Simulated  Catchments 

In  Subsection  3. 10.2  it  was  shown  that  the  mean  values  of  exterior  and  interior  link 
lengths  were  not  statistically  different.  Smart  in  his  paper  [38]  also  suggested  different 
distributions  for  areas  draining  to  an  exterior  and  interior  link.  A  t-test  similar  to  the  test 


Mean  Link  Length 
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Figure  3-50:  Link  Length  vs.  Magnitude. 
Cases  10,34,35,36 
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Figure  3-51:  Link  Length  vs.  Downstream  Magnitude. 
Cases  10,34,35,36 
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performed  in  Section  3.10.2  was  done  using  mean  contributing  areas  ( Ae  and  A;  for  exterior 
and  interior  areas  respectively)  instead  of  link  lengths.  The  values  of  Ae,  A(  are  shown  in 
Table  3-FV.  The  t-test  at  the  5%  level  of  significance  shows  AJA(  to  be  larger  than  1.  This 
result  is  in  agreement  with  [47].  The  difference  in  the  mean  of  the  area  distributions  comes 
from  the  source  area  that  contributes  only  to  exterior  links.  Given  that  no  statistical 
difference  was  found  between  mean  exterior  and  interior  link  lengths  and  that  lateral  area 
and  link  length  are  strongly  correlated,  then  the  difference  between  Ae  and  Ai  comes 
basically  from  source  areas.  The  mean  proportion  of  source  areas  to  total  area  is  27%  which 
is  a  little  low  when  compared  to  the  value  of  the  proportion  calculated  from  digital 
elevation  maps  where  its  value  is  around  35-40%  (Moglen  and  Bras,  unpublished  data,). 

3.10.4  Fundamental  Length  Scales  in  Simulated  Catchments 

In  this  subsection  commonly  studied  catchment  length  scales  will  be  measured  in  the 
resulting  simulated  networks  generated  with  the  subsurface  saturation  modification  to  the 
WBR  model.  These  measures  address  the  problem  of  determining  a  meaningful  definition 
of  mean  hillslope  length.  The  lengths  to  be  studied  are:  the  mean  link  lengths  of  exterior, 
interior  and  the  whole  set  of  links  (L€,  L,  and  L,  respectively);  the  drainage  density-based 
hillslope  length  (Lh=l/2Dd)  proposed  by  Horton  [14];  the  square  root  of  the  mean  first  order 
Strahler  area  Va,;  the  mean  value  VAfof  the  area  contributing  to  channel  heads  (called 
source  area  by  Montgomery  and  Dietrich  [28]);  and  the  lateral  hillslope  length  (defined  in 
[47]  as  the  mean  hillslope  length  of  the  lateral  contributing  areas). 

In  order  to  make  a  comparison  between  the  results  of  the  original  WBR  model  and  its 
subsurface  saturation  modification,  a  simple  linear  least-square  regression  was  fitted 
pairwise  to  the  different  length  scales  described  before.  Results  appear  in  Tables  3-V  and 
3- VI 
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VA- " 

SIMULATION 

Ae 

A; 

BR643SATU10 

15.38 

7.07 

2.34 

BR643SATU21 

8.02 

4.09 

1.91 

BR643SATU29 

20.11 

9.40 

2.26 

BR643SATU30 

10.27 

7.42 

0.75 

BR643SATU31 

9.63 

4.00 

2.80 

BR643SATU32 

14.27 

7.20 

1.79 

BR643SATU33 

10.82 

6.63 

1.25 

BR643SATU34 

18.00 

7.11 

3.04 

BR643SATU35 

16.69 

9.78 

1.41 

BR643SATU36 

12.85 

7.37 

1.47 

BR643SATU37 

18.16 

9.10 

1.99 

BR643SATU39 

13.31 

6.83 

1.89 

BR643SATU40 

51.00 

36.44 

0.79 

BR643SATU42 

2C.98 

5.86 

5.12 

BR643SATU43 

8.97 

4.45 

2.02 

BR643SATU44 

53.11 

25.06 

2.22 

BR643SATU45 

6.76 

4.50 

0.99 

Table  3-IV:  t-Test  for  Mean  Exterior  and  Interior  Link  Areas 


The  correlation  between  the  different  length  scales  is  generally  good.  The  only 
exceptions  are  the  correlation  values  between  the  mean  exterior  link  length  and  the  other 
scales,  even  though  the  correlation  between  Lt  and  VA,is  high.  In  the  WBR  model  the 
variable  with  less  correlation  against  the  others  was  the  mean  interior  link  length.  Notice 
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also  that  the  mean  first  Strahler  order  area  is  approximately  75%  larger  than  the  mean 
source  area  and  this  influences  the  relation  between  mean  and  exterior  areas  that  were 
studied  in  the  Subsection  3.10.3. 


X\Y 

L, 

Lf 

L, 

L* 

>1 

A, 

K 

L, 

1.00 

0.62 

0.96 

1.32 

0.79 

0.92 

1.99 

0.62 

1.00 

1.05 

1.82 

0.85 

1.23 

2.64 

L, 

0.96 

1.05 

1.00 

1.56 

0.82 

1.07 

2.31 

L* 

1.32 

1.82 

1.56 

1.00 

0.43 

0.61 

1.36 

>1 

0.79 

0.85 

0.82 

0.43 

1.00 

1.32 

2.84 

0.92 

1.23 

1.07 

0.61 

1.32 

1.00 

2.13 

L  ih 

1.99 

2.64 

2.31 

1.36 

2.84 

2.13 

1.00 

Table  3-V:  Linear  Regression  Between  Fundamental  Length  scales 


X\Y 

L, 

Lf 

L* 

^A, 

A, 

L 

L, 

1.00 

0.67 

0.91 

0.71 

0.90 

0.77 

0.76 

L , 

0.67 

1.00 

0.92 

0.91 

0.90 

0.96 

0.93 

Li 

0.91 

0.92 

1.00 

0.89 

0.99 

0.95 

0.93 

l* 

0.71 

0.91 

0.89 

1.00 

0.92 

0.95 

0.96 

VAi 

0.90 

0.90 

0.99 

0.92 

1.00 

0.97 

0.95 

A, 

0.77 

0.96 

0.95 

0.95 

0.97 

1.00 

0.97 

Lft 

0.76 

0.93 

0.93 

0.96 

0.95 

0.97 

1.00 

Table  3- VI:  Correlation  Coefficients  Between  Fundamental  Length  scales 
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3.11  The  Influence  of  the  Probabilistic  Distribution  of  Rainfall 

In  Section  3.2  the  excess  water  at  every  node  was  calculated  as  the  difference 
between  the  rainfall  P3  at  the  mean  flood  event  and  the  deficit  D(  calculated  using  Equation 
(3.10).  The  sum  along  the  elevation  gradient  of  the  excess  water  at  each  point  gives  the 
contributing  overland  flow  as  in  Equation  (3.14).  In  all  the  simulations  examined  up  to  this 
point,  the  calculations  of  the  amount  of  overland  flow  at  every  node  were  based  on  the 
fixed  mean  value  p3.  This  implies  that  in  the  calculations  some  nodes  were  not  saturated 
and  therefore  no  erosion  due  to  overland  flow  is  possible.  However,  given  that  the 
timescales  considered  in  the  simulations  are  very  long,  typically  on  the  order  of  millions  of 
years,  it  would  be  expected  that  some  rainfall  events  corresponding  to  floods  above  the 
mean  value  p3  would  produce  erosion  on  those  nodes  considered  as  unsaturated  before. 
Following  Willgoose  (personal  communication)  and  [19],  a  way  to  modify  the  overland 
flow  calculation  procedure  is  to  assume  a  probabilistic  distribution  over  the  rainfall 
amounts  and  calculate  the  excess  water  based  on  this  distribution. 

If  an  exponential  distribution  of  rainfall  amounts  R  with  mean  P3  is  assumed,  then  the 
excess  water  at  a  node  i  is: 

£.  =  f~  — —2^  e-^3  dR  =  i-  f"  (R-D)  e"*^  dR  =  p3  t~DJh  (3.45) 
J0  P3  P3  J£>. 

where  again  <x>  stands  for  max(x,0). 

A  node  is  saturated  if  D,  is  equal  to  0  and  in  this  case  E,=P3  as  in  the  original 
formulation.  The  difference  in  the  overland  flow  distribution  in  the  new  formulation  comes 
from  the  originally  unsaturated  nodes  which  now  have  some  excess  water  production. 
Equation  (3.45)  can  be  seen  as  a  probabilistic  proportion  of  overland  flow  assigned  to  the 
nodes.  Notice  however,  that  in  this  approach,  even  though  the  probabilistic  distribution  of  R 
is  taken  into  account,  no  stochastic  term  is  included  in  the  model’s  formulation. 
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Furthermore,  because  of  the  special  form  of  the  exponential  distribution,  the  calculation  of 
E(  still  depends  on  the  parameters  of  the  original  model. 

The  overland  flow  at  a  node  i  is  now: 

Q,  -£,■  + 1 Q,  >>,  (3-46) 

1 

and  the  other  equations  of  the  model  remain  the  same. 

Using  the  same  parameters  and  the  initial  elevation  field  of  the  simulation  example 
presented  in  Section  3.4,  a  new  simulation  including  the  approach  of  this  section  was 
performed.  Comparisons  between  both  runs  will  now  be  presented.  Figure  3-52  shows  the 
final  networks  for  the  original  simulation  BR643SATU10  and  the  modified  version 
BR661SATU10.  The  differences  are  small  but  noticeable.  Figure  3-53  to  3-56  show  the 
evolution  of  the  elevation  distribution  with  time.  Figures  3-57  and  3-58  correspond  to  the 
evolution  in  time  of  the  drainage  density  and  magnitude  of  the  network. 

Finally,  Figure  3-59  shows  the  evolution  of  the  hypsometric  curves  for  both  cases. 
This  figure  shows  clearly  the  effect  of  the  new  overland  flow  distribution  on  the  catchment 
evolution.  It  can  be  seen  that  those  nodes  at  higher  elevation,  which  usually  fall  into  the 
group  of  unsaturated  nodes  because  of  their  smaller  areas  and  larger  slopes,  are  eroded 
faster  in  the  new  formulation.  The  sharpness  of  the  knee  is  reduced  because  of  the  erosion 
in  the  unsaturated  nodes  and  the  hypsometric  curves  in  both  cases  are  almost  identical  in 
the  zones  of  channel  and  saturated  nodes.  The  overland  flow  in  the  unsaturated  nodes,  even 
though  small,  is  enough  to  produce  changes  in  the  network.  Similar  behavior  was  observed 
in  different  simulations  that  used  the  same  parameter  values  as  the  experiments  described  in 
previous  sections.  Figure  3-60  shows  the  hypsometric  curve  evolution  for  case  33  that  was 
analyzed  in  Section  3.6. 

At  the  beginning  of  the  basin  evolution  almost  all  the  nodes  are  saturated  and  no 
difference  exists  between  both  simulations.  As  time  passes,  some  nodes  become 
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Figure  3-56:  Elevations  at  Time  7000.  Case  BR661SATU10 


Figure  3-57:  Drainage  Density  Evolution.  Cases  BR643SATU10  and  BR661SATU10 


-87- 


Figure  3-58:  Magnitude  Evolution.  Cases  BR643SATU10  and  BR661SATU10 


Figure  3-59:  Hypsometric  Curve  Evolution.  Cases  BR643SATU10  and  BR661SATU10 
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Figure  3-60:  Hypsometric  Curve  Evolution.  Cases  BR643SATU33  and  BR661SATU33 

unsaturated  and  do  not  contribute  to  overland  flow  in  the  original  formulation.  However, 
even  though  excess  water  is  produced  in  those  nodes  in  the  second  fomulation,  the 
exponential  form  of  the  distribution  makes  this  quantity  small.  For  example,  in  simulation 
BR661SATU10  the  probabilistic  proportion  of  overland  How  in  unsaturated  nodes  was  one 
to  three  orders  of  magnitude  smaller  than  in  saturated  nodes.  The  effect  of  E,  is  small  in  the 
overall  evolution  because  in  both  cases  the  evolution  of  the  catchments  is  led  by  the 
channelized  nodes. 

Table  3- VII  compares  the  statistics  of  the  network  for  the  cases  BR643SATL10  and 
BR661SATU10. 

The  role  of  T.  also  changes,  being  more  important  in  the  new  approach  described  in 
this  section.  In  the  original  formulation  T.,  through  its  effect  on  D  ,  determined  which  nodes 
were  saturated  and  which  were  not.  T.  was  also  involved  in  the  overland  How  calculation  in 
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Statistic 

BR643SATU10 

BR661SATU10 

RfcU-2) 

5.08 

4.62 

Rs(l-2) 

2.85 

2.69 

Rl(1-2) 

2.60 

2.91 

R*(l-2) 

5.85 

6.14 

R*(2-3) 

4.00 

4.33 

Rs(2-3) 

2.15 

2.02 

Rl(2-3) 

1.82 

1.64 

R„(2-3) 

4.96 

4.63 

R* 

4.75 

4.47 

Rs 

2.67 

2.53 

R  L 

2.35 

2.57 

R, 

5.53 

5.73 

K 

1.96 

1.68 

ei 

2.67 

2.26 

d; 

11.54 

11.73 

Magnitude 

61 

60 

Mean  Relief 

9.90 

9.90 

Table  3-VII:  Final  Network  Statistic  Values. 

Simulations  BR643SATU10  and  BR661SATU10 

those  nodes  that  were  near  complete  saturation  (i.e.  those  nodes  with  D>0  but  P3-D>0).  In 
the  approach  described  in  this  section,  Tc  is  important  in  calculating  the  overland  flow  at 
every  unsaturated  node  through  Equation  (3.45V  Tc  affects  the  probabilistic  proportion  of 
overland  flow  assigned  to  the  unsaturated  nodes.  Figures  3-61  to  3-65  show  simulations 
with  increasing  values  of  TSC.  The  elevation  fields  and  the  resulting  networks  are 
presented.  The  effect  of  Tc  can  be  clearly  seen  in  the  variation  of  the  hypsometric  curves. 
Figures  3-66  to  3-68  correspond  to  simulations  with  TSC  equal  to  0.276xl03,  0.920xl03 
and  4.600xl03.  Notice  how  the  knee  in  the  hypsometric  curves  tends  to  reduce  as  TSC 
decreases.  The  increasing  number  of  saturated  nodes  and  the  increasing  overland  flow 
production  in  unsaturated  nodes  are  the  reasons  for  this  behavior.  In  the  limit,  as  TSC 
decreases,  the  subsurface  saturation  modification  tends  to  the  original  WBR  model. 
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Figure  3-62:  TSC=0.552xl03 
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Figure  3-67:  Hypsometric  Curve  Evolution  TSC=0.920xl03 


Figure  3-68:  Hypsometric  Curve  Evolution  TSC=4.60xl03 
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Chapter  4 

Sensitivity  to  Initial  Conditions  in  the  WBR  model 

4.1  Introduction 

The  networks  produced  b\  both  the  WBR  model  and  its  subsurface  saturation 
modification  have  a  random  appearance.  Keeping  all  the  parameter  values  constant,  and 
varying  the  initial  elevation  field  by  a  small  amount  produce  completely  different  networks. 
Analytical  results  on  the  Meinhardt  equations  ''described  in  Section  2.3.)  show  that  small 
perturbations  in  the  system  can  grow  and  produce  completely  different  r.  ^ults,  making  the 
system  chaotic;  it  is  this  transient  chaos  that  causes  the  random  appearance.  The  special 
fonm  of  the  Meinhardt  system  allows  an  analytical  derivation  of  these  facts,  but  this 
procedure  is  not  directly  applicable  to  the  WBR  catchment  evolution  model.  The 
integration  of  contributing  How  to  each  nodt  along  the  elevation  gradient  is  one  of  the 
major  difficulties.  Nevertheless,  numerical  exneriments  can  be  performed  in  order  to 
establish  whether  the  chaotic  behavior  found  in  the  Meinhardt  system  is  similar  to  the  WBR 
model  and  to  what  extent. 

The  organization  of  the  chapter  is  as  follows:  a  short  review  of  definitions  and 
common  procedures  in  the  study  of  chaotic  dynamical  systems  is  presented.  Then,  after 
defining  an  appropiate  measure  between  elevation  fields  called  the  relief  difference,  various 
experiments  are  presented.  All  these  experiments  try  to  examine  the  validity  of  the 
hypothesis  that  evolutions  that  begin  at  very  similar  initial  conditions  evolve  in  a 
completely  different  manner.  Finally,  the  influence  of  the  different  parameters  of  the  WBR 
model  over  the  dynamical  evolution  of  the  system  is  examined. 
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4.2  Sensitivity  to  Initial  Conditions  and  Lyapunov  Exponents 

There  is  still  no  agreement  for  an  exact  definition  of  a  chaotic  dynamical  system. 
Further  research  and  stronger  mathematical  theory  are  required  before  a  consensus  about 
such  definition  is  reached.  However,  there  is  one  property  that  all  chaotic  systems  share  and 
which  makes  them  interesting.  This  property  is  the  sensitivity  to  initial  conditions 
(henceforth  called  SIC).  The  SIC  property  means  that  nearby  trajectories  in  the  phase  space 
will  move  away  exponentially  fast  from  each  other  on  the  average.  It  implies  then  that  no 
matter  how  precisely  variable  measurements  are  taken  in  a  chaotic  system,  errors  will  grow 
and  dominate  the  solution  making  all  long-term  predictability  impossible. 

Complexity  is  by  no  means  a  requirement  in  order  for  a  system  to  exhibit  SIC.  A 
system  as  simple  as  the  logistic  equation  xfl+|=axn(l-xn)  (with  a  a  parameter  in  [0,4]  and  \n 
in  [0,1])  can  exhibit  this  property.  Figure  4-1  shows  the  evolution  of  two  initial  conditions 
x0=0.2  and  x0=0.2001  in  different  regimes:  one  periodic  and  one  chaotic.  In  the  periodic 
regime  (a=3.5)  both  solutions  evolve  side  by  side  and  the  distance  between  them  remains 
bounded  for  all  time.  On  the  other  hand,  in  the  chaotic  regime  (a=4.0)  even  though  the 
difference  between  the  two  initial  conditions  is  very  small,  after  some  time  the  solutions 
evolve  in  a  completely  different  fashion.  Any  small  error  can  produce  loss  of  predictability 
in  the  long  term. 

The  rates  at  which  nearby  solutions  diverge  are  called  Lyapunov  exponents.  Their 
sign  indicates  whether  nearby  trajectories  diverge  exponentially  (i.e.  if  SIC  is  present  or 
not);  positive  signs  imply  chaotic  behavior.  Lyapunov  exponents  also  indicate  the  time 
scale  on  which  predictions  within  a  certain  accuracy  level  can  be  performed. 

In  mathematical  terms  [27],  if  d0  is  a  measure  of  the  separation  between  two 
trajectories  of  a  dynamical  system  at  a  certain  time,  the  distance  between  them  at  a  later 
time  can  be  described  as: 


'  '  I  - - 1 - - - 1 - - - 1 - -  I 
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Tune  n 


Figure  4*1:  Logistic  Equation.  a=3.5  and  4.0  respectively 
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d(t)  =  dQ  eh  (4-D 

where  X,  is  the  Lyapunov  exponent  Figure  4-2  shows  this  property  graphically.  If  X  is 
positive,  the  system  displays  SIC  and  behaves  chaotically,  otherwise  the  system  behaves 
regularly.  However,  it  is  required  that  the  dynamical  system  has  dimensionality  three  or 
larger  in  order  to  exist  the  possibility  of  chaotic  behavior. 


Figure  4-2:  Exponential  Separation  of  Trajectories  in  a  Chaotic  System 

Of  course.  Equation  (4.1)  is  not  valid  for  all  time  because  the  system  is  bounded  and 
the  distance  between  trajectories  cannot  go  to  infinity.  This  is  the  reason  why  some  average 
is  required.  Following  Wolf  et  al.  [48],  a  way  to  calculate  the  largest  Lyapunov  exponent  is 
to  take  two  nearby  initial  conditions,  follow  their  trajectories  in  time  and  when  the  distance 
between  them  becomes  too  large,  to  replace  one  of  the  points  with  another  one  near  the 
base  trajectory.  Section  4.6  will  examine  this  idea  in  detail. 


The  Lyapunov  exponent  is  calculated  then  as: 

,  1  i  .  <*('*) 

X  =  —  y  log— — - 

^0  ks  1 


(4.2) 


where  t(  is  the  time  at  which  a  new  measurement  is  taken,  d(tt)  the  distance  between 
trajectories  at  time  t*  and  d^t*)  the  initial  distance  between  the  base  trajectory  and  the  new 
trajectory  chosen  at  time  t*. 


If  the  Lyapunov  exponent  X.  is  going  to  describe  long-term  behavior  and  have 
dependence  on  the  initial  condition  used  for  its  calculation,  the  average  should  continue  for 
a  long  enough  time.  Notice  also,  if  the  system  is  1-D  and  the  time  is  small,  the 
relation  tends  to  the  derivative /'(**_,)  of  the  solution  of  the  system. 

For  systems  of  dimensionality  larger  than  one,  a  whole  set  of  Lyapunov  exponents 
can  be  defined,  each  of  which  describes  the  divergence  of  trajectories  in  a  particular 
direction.  In  a  n-dimensional  system  a  small  n-ball  of  initial  conditions  evolves  into  an  n- 
ellipsoid  with  time.  The  exponential  rate  at  which  the  different  n  principal  axis  of  this 
ellipsoid  grows  or  decays  is  the  generalization  of  the  concept  of  the  1-D  Lyapunov 
exponent.  Figure  4-3  from  [27]  explains  this  effect.  The  sum  of  all  the  Lyapunov  exponents 
is  the  time-averaged  divergence  of  its  phase  space.  Wolf  et  ai.  [48]  use  this  idea  to  calculate 
Lyapunov  exponents:  the  growth  of  length  elements  (distance  between  two  trajectories) 
gives  the  largest  Lyapunov  exponent;  the  growth  of  area  elements  (using  three  trajectories) 
gives  the  sum  of  the  largest  two  exponents  and  so  on. 

Tim* 

Volume 

Figure  4-3:  The  Divergence  of  Volumes  in  Chaotic  Systems 

A  word  of  caution  is  necessary  at  this  point:  there  is  a  stricter  definition  of  Lyapunov 
exponents  and  expansion  rates  of  a  dynamical  system.  This  definition  examines  the 
eigenvalues  of  the  transformation  matrix  corresponding  to  the  linearization  of  the  original 
dynamical  system.  Given  that  the  linearization  as  well  as  the  direction  of  its  eigenvectors 
vary  with  time,  an  average  over  time  of  the  eigenvalues  is  used  to  define  Lyapunov 
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exponents  (see  Berge  et  al.  [3]  for  details).  In  this  case  the  exponents  also  measure  average 
linearized  expansions  and  contractions.  It  is  expected  that  both  definitions  (the  experimental 
and  the  theoretical)  are  strongly  related  but  further  research  is  still  necessary. 

4.3  SIC  in  the  Catchment  Evolution  Model 


The  WBR  model  has  not  only  non-linear  time  evolution  but  it  is  also  a  spatially 
extended  system.  Apart  from  some  works  by  Kaneko  and  co-workers  [16],  [17]  and  [23]  in 
coupled  logistic  maps,  studies  of  chaotic  behavior  in  spatially  extended  systems  do  not 
exist  They  are  to  be  expected  in  future  years  once  the  mathematical  theory  advances 
enough  to  look  at  problems  of  turbulence  where  both  time  and  space  are  inherently 
unpredictable.  This  section  will  try  to  apply  and  generalize  the  concepts  described  in 
Section  4.2  to  the  WBR  model. 

The  first  problem  is  to  define  a  measure  of  distance  between  trajectories  in  order  to 
study  the  evolution  of  d(t)  as  in  Equation  (4.1).  A  plausible  definition  of  the  distance 
between  two  elevation  fields  Z,  and  Z2  is: 

dn(t)  *  IIZj(f)-Z2(f)ll  =  £  (4.3) 

allnodu 

where  is  the  elevation  of  node  (ij)  at  time  t.  Measure  (4.3)  will  be  called  the  relief 
difference  between  elevation  fields  Z,  and  Z2.  In  stricter  mathematical  terms,  if  the 
elevation  at  each  node  is  considered  as  an  independent  component  of  an  n*m-dimensional 
vector  (where  m  and  n  are  the  vertical  and  horizontal  dimensions  of  the  numerical  grid 
respectively),  then  (4.3)  is  the  L0-noim  in  the  n*m-D  space.  Of  course,  any  other  valid  norm 
in  this  space,  for  example  the  square  root  of  the  sum  of  the  squares  of  the  elevation 
differences  at  grid  nodes  can  also  be  used. 

The  SIC  property  can  be  tested  then,  under  different  conditions  in  the  WBR  model. 
Two  initial  elevation  fields  are  taken  such  that  the  distance  (i.e.  the  relief  difference) 
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between  them  is  small  and  the  evolution  of  both  catchments  is  followed  measuring  dX2(t).  If 
the  system  were  regular  d12(r)  would  remain  bounded  for  ail  times.  In  the  experiments 
performed,  a  base  elevation  field  Zb  was  chosen  randomly:  elevations  were  drawn  from  an 
uniform  distribution  with  a  very  small  coefficient  of  variation.  Its  evolution  was  compared 
against  trajectories  of  points  in  phase  space  (initial  elevation  fields  in  real  space)  near  the 
base  elevation  field.  The  details  of  how  these  points  are  chosen  is  described  in  the  following 
two  sections. 

4.4  Behavior  of  One-Node  Perturbations 

Once  the  base  initial  elevation  field  Zj(0)  is  fixed,  a  way  to  choose  another  elevation 
field  Z2  "near"  Zx  is  to  perturb  the  elevation  at  one  of  the  nodes  by  a  small  amount.  The 
measure  of  how  separate  Z2  is  from  Zx  is  given  by  the  relief  difference  between  them. 
Given  that  only  the  elevation  at  one  node  is  varied,  d,2(0)  is  equal  to  th  variation  in  z  at 
that  node.  Figure  4-4  shows  the  location  of  points  at  which  the  elevation  was  perturbed  in 
each  case.  The  final  network  is  presented  only  for  illustrative  purposes;  at  the  beginning  of 
the  simulation  no  channel  network  exists.  The  elevation  was  changed  by  the  same  amount 
in  every  experiment  (approximately  3%  of  the  notch).  Figure  4-5  shows  the  evolution  in 
time  of  the  value  of  the  relief  difference  (i.e.  the  distance  in  phase  space)  between  Zx  and 
Zj.  Notice  how,  for  some  of  the  cases,  d12(t)  grows  exponentially  with  time  and  continue 
this  expansion  until  the  network  captures  the  entire  catchment.  Once  the  catchment  is 
captured,  at  approximately  dimensionless  time  1500  in  the  experiment  shown,  the  growth  of 
d12(t)  stops.  Then,  the  relief  difference  between  the  random  fields  begins  to  decrease  slowly 
because  in  the  case  with  episodic  uplift  both  catchments  tend  ultimately  to  a  flat  plane.  The 
exponential  separation  between  trajectories  is  an  indication  of  the  transient  chaos  present  in 
the  WBR  model.  The  logarithmic  scale  in  the  relief  difference  axis  makes  it  possible  to 
infer  visually  that  the  sign  of  the  Lyapunov  exponent  (which  corresponds  to  the  slope  of  the 
graphs)  is  positive. 


Figure  4-4:  Location  of  Nodes  with  Perturbed  Elevations 


Figure  4-5:  Difference  Relief  for  One-Node  Perturbed  Elevation  Fields 
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Some  other  trajectories  do  not  diverge  exponentially  from  the  base  trajectory  during 
the  whole  network  expansion  period  as  can  be  seen  in  Figure  4-5.  In  these  cases  d12(t) 
decreases  rapidly  even  if  there  was  a  small  initial  expansion.  The  interesting  fact  is  that  the 
final  networks  are  identical  in  the  bounded  case  while  on  the  exponentially  divergent  cases 
the  resulting  networks  are  completely  different.  In  terms  of  the  phase  space,  the  base 
trajectory  attracts  some  of  the  perturbed  trajectories.  The  relationship  between  exponential 
divergence  and  different  final  network  is  one-to-one.  Final  networks  are  presented  in  Figure 
4-6 

There  are  some  features  common  to  the  chaotic  behavior  of  the  Meinhardt  and  the 
WBR  systems.  The  analytical  formulation  of  the  Meinhardt  model  permits  the  ellucidation 
of  the  way  in  which  perturbations  grow  and  dominate  the  solution  (see  [47]).  Fluctuations 
in  the  channel  formation  function  near  the  branch  will  propagate  in  an  unstable  manner,  and 
fluctuations  away  from  the  differentiation  process  are  reduced.  However,  if  the 
differentiation  boundary  passes  at  a  later  time  through  the  zone  where  traces  of  the  original 
fluctuations  still  exist,  then  they  also  propagate  unstably.  In  the  WBR  model,  elevation 
perturbations  affect  both  flow  directions  and  channel  formation  function  values  (the  latter 
through  the  slope  factor  in  the  channel  formation  function  equation).These  variations  enter 
into  the  model  equations  and  eventually  dominate  the  solution.  Notice  also  that 
perturbations  at  nodes  away  from  the  differentiation  process  do  not  grow  at  all. 

This  unstable  propagation  of  fluctuations  raises  the  issue  of  whether  the  numerical 
solution  of  the  WBR  model  is  really  a  true  catchment  evolution.  Given  that  numerical 
errors  grow  unstably,  the  numerical  solution  may  not  be  a  true  trajectory.  A  possibility  is 
that  the  numerical  solution  bounces  around  different  real  trajectories  of  the  dynamical 
system.  However,  recent  research  by  Yorke,  Grebogi  and  coworkers  [12]  in  simple  systems 
seems  to  indicate  that  there  is  always  a  true  solution  near  a  numerical  one  in  chaotic 
dynamical  systems.  Even  though  further  research  is  still  necessary,  these  preliminary  steps 


Figure  4-6:  Final  Networks  for  One-Node  Perturbed  Elevation  Fields 
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give  confidence  on  the  good  behavior  of  numerical  solutions  of  chaotic  dynamical 
solutions.  Of  course,  long-term  predictability  is  not  possible  in  these  systems  if  the 
precision  is  finite. 

4.5  Expansion  of  a  Ball  in  the  Phase  Space  of  Elevation  Fields 

In  Section  4.4  Z,  was  constructed  by  changing  the  elevation  at  one  node  in  Z,. 
Another  possibility  is  to  change  the  elevations  at  every  node,  i.e.  to  take  a  completely 
different  initial  elevation  field  drawn  from  the  same  distribution  to  which  Z,  belongs.  Given 
that  the  coefficient  of  variation  of  the  uniform  distribution  of  node  elevations  is  small,  the 
relief  difference  d12(0)  will  be  also  small.  Thinking  in  terms  of  the  n*m-D  phase  space, 
perturbations  at  one  node  correspond  to  movement  along  one  of  the  axes  while 
perturbations  corresponding  to  new  random  elevation  fields  are  represented  by  points 
positioned  away  from  the  origin  into  the  space.  Various  trajectories  corresponding  to 
different  initial  random  fields  were  compared  in  time  with  the  evolution  of  Zb.  These  initial 
elevation  fields  can  be  viewed  as  a  ball  in  the  phase  space.  The  deformation  of  this  small 
ball  was  followed.  The  goal  was  to  examine  the  possibility  of  exponential  growth  of  the  ball 
as  in  Figure  4-3  in  section  4.2  .  If  exponential  growth  occurs,  the  SIC  property  is  shown  not 
to  be  restricted  to  single  node  perturbations  and  the  sum  of  the  Lyapunov  exponents  can  be 
studied. 

Figure  4-7  shows  three  measures  of  the  ball’s  evolution:  the  maximum  and  the 
minimum  values  of  the  set  of  the  relief  differences,  and  the  maximum  distance  between 
trajectories.  All  three  measures  provide  insight  into  the  ball  expansion  process.  The  second 
measure  proves  that  it  is  not  a  small  group  of  trajectories  that  is  diverging  from  Z,  but  that 
effectively  all  of  them  are  moving  away  exponentially  fast.  The  last  measure  gives  an  idea 
of  the  largest  principal  axis  size  (its  exponential  rate  is  equivalent  to  the  largest  Lyapunov 
exponent,  see  [27]). 
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Figure  4-7:  Three  Measures  of  the  Expansion  of  a  Volume  of 
Nearby  Initial  Random  Fields  in  Phase  Space 

Even  though  a  very  large  number  of  points  is  required  to  determine  with  precision  the 
value  of  the  sum  of  the  Lyapunov  exponents  based  on  the  volume  expansion  of  the  ball  of 
initial  conditions.  Figure  4-7  gives  adequate  certainty  about  the  positive  sign  of  the  largest 
Lyapunov  exponent  and  therefore  the  chaotic  nature  of  the  WBR  system. 


4.6  Another  Test  Related  to  the  SIC  Hypothesis 


Another  test,  based  on  ideas  by  Wolf  et  al.  (48]  is  to  examine  not  only  the  evolution 
of  trajectories  corresponding  to  initial  points  near  the  base  elevation  field,  but  also  to  look 
at  evolutions  that  begin  near  the  base  trajectory  at  later  times,  as  was  shown  in  Figure  4-2. 
This  method  has  a  clear  advantage:  given  that  the  size  of  the  domain  under  consideration  in 
the  phase  space  is  finite,  then  exponential  divergence  cannot  continue  forever  (as  can  be 
seen  in  Figure  4-5).  That  is  the  reason  why  in  chaotic  dynamical  systems  both  time  and 
space  averages  have  to  be  taken  in  order  to  calculate  Lyapunov  exponents. 
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For  this  reason  perturbations  of  the  base  elevation  field  were  performed  at  different 
times  and  their  trajectories  were  compared  to  the  base  field  evolution.  Results  appear  in 
Figure  4-8.  In  this  figure  the  perturbations  that  begin  at  different  rimesteps  have  been 
translated  to  coincide  with  the  origin.  Times  in  the  legend  indicate  the  step  at  which  the 
elevation  field  is  perturbed.  Notice  that  SIC  and  chaotic  behavior  is  present  only  while  the 
network  is  still  growing.  The  network  grows  up  to  nondimensionai  time  2000;  perturbations 
initiated  after  this  point  decay  in  time. 

Once  the  network  captures  the  catchment  and  provides  rigidity  to  the  elevation  field, 
perturbations  decrease  with  time.  The  WBR  system  exhibits  what  is  called  transient  chaos. 
The  largest  Lyapunov  exponent  can  be  calculated  using  the  different  perturbations  and 
Equation  (4.2). 


Figure  4-8:  Relief  Difference  of  Nearby  Points 
Along  the  Base  Trajectory 
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4.7  One>Node  Perturbation  Thresholds 

By  performing  various  experiments  like  the  ones  described  in  Section  4.4  it  was 
found  that  there  was  a  perturbation  size  threshold  that  has  to  be  exceeded  in  order  to  have 
exponential  divergence  of  trajectories  in  the  system.  This  is  to  due  to  the  one-to-one 
relation  between  the  exponential  divergence  of  trajectories  and  the  change  of  network  form, 
as  well  as  the  existence  of  an  channel  formation  function  threshold  that  has  to  be  surpassed 
in  order  to  trigger  channelization.  Of  course,  the  fact  that  only  eight  possible  flow  directions 
are  considered  at  each  node  also  influences  the  perturbation  threshold. 

In  classical  chaotic  studies  the  SIC  property  is  always  thought  of  in  terms  of 
infinitesimal  perturbations.  The  WBR  model  behaves  differently  on  this  respect.  Consider 
the  base  elevation  field  Zb  and  choose  a  node.  If  perturbations  are  made  to  the  elevation  at 
that  node,  the  perturbed  fields  are  located  along  one  of  the  axes  of  the  phase  space.  The 
evolution  of  the  relief  difference  for  perturbations  performed  this  way  over  Zb  at  different 
nodes  are  shown  in  figures  4-9  and  4-10. 

Notice  that  if  the  perturbation  is  small  enough,  even  though  the  separation  grows 
exponentially  for  some  time,  the  distance  between  trajectories  stops  growing;  it  decreases 
even  before  the  network  finishes  the  capture  process.  There  exists  a  threshold  value, 
dependent  on  the  elevation  at  neighboring  nodes  and  on  the  parameter  values  but  it  is 
certainly  small.  Above  this  threshold  the  perturbed  trajectory  continues  its  exponential 
separation.  These  are  the  two  different  behaviors  observed  in  Figures  4-9  and  4-10.  In  the 
experiments  the  original  node  elevation  was  10.10428.  if  the  elevation  at  that  node  in  Zj 
were  larger  than  10.10850  or  smaller  than  10.10258  then  the  exponential  separation  occurs 
during  the  whole  network  expansion  period.  The  jump  at  these  two  values  occurs,  not 
surprisingly  with  the  change  in  network  form.  Regions  in  the  catchment  begin  to  evolve  in  a 
different  manner  and  the  trajectories  in  phase  space  diverge  from  each  other.  Thinking  in 
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terms  of  the  phase  space,  a  cylinder  of  very  small  size  exists  around  the  base  trajectory  in 
which  nearby  trajectories  behave  differently  from  those  outside  the  cylinder. 


Figure  4-9:  Perturbation  Threshold  One-Node  Perturbations. 
Upward  Perturbations 
Original  Node  Elevation  10.10428 


The  figures  referred  to  earlier  show  a  discontinuous  jump  between  two  regimes  that 
vary  continuously  with  the  size  of  the  elevation  perturbation.  Figure  4- 1 1  coresponds  to  the 
slopes  of  least-squares  lines  fitted  to  the  portion  of  the  curves  in  the  time  interval 
[300,1000]  in  the  semilog  plane.  The  results  shown  in  Figure  4-11  are  based  on  both 
Figures  4-9  and  4-10.  It  can  be  observed  that  near  the  original  node  elevation  (i.e.  inside  the 
cylinder  in  phase  space)  slopes  are  negative  and  the  base  trajectory  attracts  neighbor 
trajectories. 
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Figure  4-10:  Perturbation  Threshold  One-Node  Perturbations. 
Downward  Perturbations 
Original  Node  Elevation  10.10428 


Figure  4-11:  Slopes  of  Relief  Difference  Curves 
Original  Node  Elevation  10.10428 
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4.8  Parameter  Dependence  of  SIC 

A  number  of  parameters  can  be  varied  in  the  WBR  model,  and  when  this  occurs  the 
chaotic  behavior  of  the  system  changes  as  well.  The  most  important  parameters  to  be 
studied  are  the  exponents  of  the  flow  and  slope  terms  in  the  sediment  transport  and  channel 
formation  function  equations  (m^n^nij,  and  n5  respectively),  the  channel  formation 
function  threshold  (l/cj),  the  coefficient  (33  used  in  the  relation  flow-area  (a  rainfall 
measure),  the  coefficient  (Sj  in  the  sediment  transport  formula  and  the  relation  between 
channel  and  hillslope  sediment  transport  (0(). 

In  the  next  set  of  experiments  the  evolutions  of  a  base  elevation  field  Zb  and  a 
perturbed  one  Z2  are  again  followed.  The  initial  elevation  field  Z2(0)  was  chosen  as  in 
Section  4.5.  This  time,  however,  the  fields  under  observation  are  kept  constant  and  the 
parameter  values  are  varied.  This  procedure  will  show  the  dependence  of  the  divergence  of 
the  trajectories  with  respect  to  the  parameter  values.  The  following  graphs  show  the 
exponential  trajectory  divergence  of  the  same  two  random  initial  elevation  fields  for 
different  parameter  values.  In  each  case  a  semilog  graph  corresponding  to  the  long-term 
behavior,  a  semilog  graph  that  covers  the  network  growth  period,  and  an  arithmetic  graph 
coresponding  to  that  same  period  are  presented.  The  first  graph  illustrates  the  behavior  of 
the  system  during  the  dying  period  of  the  catchment,  the  second  graph  shows  the  whole 
period  of  network  growth  and  the  third  graph  presents  in  detail  the  most  intense  period  of 
network  expansion. 

An  increasing  trend  of  the  exponential  rate  of  separation  with  increasing  values  of  m, 
is  observed  in  Figures  4-12  to  4-14.  This  effect  appears  simultaneously  with  an  increase  in 
drainage  density  in  the  final  networks:  the  perturbation  information  is  communicated  over 
the  whole  catchment  by  the  network.  Without  network  there  is  no  spatial  communication  of 
elevation  differences  and  no  exponential  divergence  exists.  The  same  effect  is  observed  in 
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the  channel  formation  function  threshold  case  (Figures  4-15  to  4-17).  As  c,  increases,  the 
drainage  density  increases  because  the  channel  formation  function  threshold  is  reduced  and 
the  exponential  rate  of  divergence  increases.  A  change  of  behavior  occurs  around  Cj=0.0007 
in  this  experiment  because  the  network  begins  to  lose  its  branched  structure  and  a  "blob"  is 
formed  instead. 


Figure  4-12:  Relief  Difference  Under  mj  Variations 

The  variation  with  n1  is  in  the  opposite  direction.  In  Figures  4-18  to  4-20,  as  n, 
increases,  sharp  peaks  and  steep  slopes  tend  to  be  flattened  out  faster  because  the  sediment 
transport  is  larger  everything  else  being  the  same,  so  the  separation  rate  between  trajectories 
is  smaller.  The  variations  with  P3  shown  in  Figures  4-33  to  4-35  correspond  to  the  relation 
between  Lyapunov  exponents  and  rainfall.  In  this  case  the  drainage  density  is  also  a  key 
factor  in  explaining  the  behavior  observed.  Figures  4-24  to  4-26  correspond  to  O, 
variations.  Notice  how  the  relief  difference  decreases  at  long  time  in  a  manner  directly 
correlated  to  the  size  of  Or 
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Figure  4-13:  Relief  Difference  Under  nij  Variations  During  Network  Growth 
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Figure  4-14:  Relief  Difference  Under  raj  Variations. 
Arithmetic  Scale 
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Figure  4-15:  Relief  Difference  Under  ct  Variations 
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Figure  4-16:  Relief  Difference  Under  cs  Variations  During  Network  Growth 
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Figure  4-17:  Relief  Difference  Under  c,  Variations. 
Arithmetic  Scale 


Figure  4-18:  Relief  Difference  Under  n,  Variations 
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Figure  4-19:  Relief  Difference  Under  nj  Variations  During  Network  Growth 


Figure  4-20:  Relief  Difference  Under  Variations. 
Arithmetic  Scale 


Figure  4-22:  Relief  Difference  Under  (53  Variations  During  Network  Growth 


Difference 


Difference  Difference 


-118- 


Figure  4-25:  Relief  Difference  Under  O,  Variations  During  Network  Growth 


Figure  4-26:  Relief  Difference  Under  O,  Variations. 
Arithmetic  Scale 
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Figure  4-30:  Relief  Difference  Under  n5  Variations 


Difference 


0 


Time  10 3 


Figure  4-31:  Relief  Difference  Under  n5  Variations  During  Network  Growth 
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Figure  4-33:  Relief  Difference  Under  P,  Variations 


Figure  4-34:  Relief  Difference  Under  3,  Variations  During  Network  Growth 
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Figure  4-35:  Relief  Difference  Under  (5,  Variations. 

Arithmetic  Scale 

Dependence  on  m5  (Figures  4-27  to  4-29)  as  well  as  on  n5  (Figures  4-30  to  4-32)  is 
very  similar  to  the  dependence  on  and  n,  respectively.  Finally,  the  variation  in  the 
multiplicative  factor  (3,  of  the  sediment  transport  equation  is  shown  in  Figures  4-33  to  4-35. 
Of  course,  some  of  the  parameter  values  may  be  out  of  the  physically  valid  range  but  this 
analysis  provides  additional  insight  on  the  behavior  of  the  model  as  a  dynamical  system  in 
itself. 

4.9  Final  Remarks 

The  sensitivity  to  initial  conditions  and  the  chaotic  behavior  of  systems  are  properties 
that  are  going  to  be  studied  in  much  more  detail  in  future  years  because  of  their 
implications  in  predictability.  In  this  chapter  it  has  been  shown  how  this  property  produces 
completely  different  networks  out  of  processes  with  identical  parameter  values  but  with 
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small  differences  in  the  initial  elevation  fields.  The  intrinsic  temporal  and  spatial  variability 
of  the  resuits  imply  the  necessity  of  a  statistical  analysis  of  the  networks.  Furthermore,  the 
implications  of  the  SIC  property  on  long-term  predictability  have  to  be  taken  into  account  if 
the  WBR  model  is  to  be  used  in  future  years  for  field  predictions. 

The  Lyapunov  exponents  were  shown  to  have  a  positive  sign  in  the  system  using 
different  tests  that  have  been  developed  in  the  literature  under  a  variety  of  conditions.  The 
fact  that  the  signs  are  positive  proves  the  existence  of  transient  chaos  in  the  evolution  of  the 
system.  However,  calculations  of  the  exact  values  of  the  Lyapunov  exponents  require  a 
very  large  number  of  experiments  or  a  stronger  theory  on  spatially  extended  systems  which 
has  not  been  developed  yet.  The  exact  values  of  these  exponents  will  provide  a  bound  for 
long-term  predictions  using  this  model  and  this  may  be  one  of  the  most  important 
applications  of  the  understanding  of  the  chaotic  behavior  of  the  system. 
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Chapter  5 
Conclusions 


5.1  Summary  of  Results 

This  work  consists  of  two  sections.  Both  sections  study  different  aspects  of  the 
channel  network  and  catchment  evolution  model  developed  by  Willgoose,  Bras  and 
Rodriguez-Iturbe  [47],  henceforth  called  the  WBR  model.  This  model  uses  a  sediment 
transport  continuity  equation  to  model  elevations,  and  concepts  of  erosion  engineering  to 
define  a  channel  formation  function  that  triggers  the  growth  of  channels.  The  overland  flow 
production  mechanism  used  in  the  WBR  model  is  the  Hortonian.  One  of  the  most 
interesting  features  of  the  WBR  model  is  the  way  hillslopes  and  channels  interact  through 
the  different  sediment  transport  processes.  The  purpose  of  the  first  section  of  this  work  is  to 
study  the  effect  of  the  overland  flow  production  mechanism  at  hillslopes  on  the  overall 
behavior  of  the  system.  This  first  pan  changes  the  Hortonian  runoff  of  the  WBR  model  to 
the  subsurface  saturation  mechanism.  In  the  Hortonian  case  overland  flow  is  assumed  to  be 
generated  uniformly  over  the  whole  area  whereas  in  the  subsurface  saturation  mechanism 
only  saturated  areas  contribute  to  the  overland  flow. 

After  a  review  of  previous  work,  a  modified  version  of  the  WBR  model,  including 
the  subsurface  saturation  runoff  mechanism,  is  presented  in  Chapter  3.  A  series  of 
simulations  were  shown  and  their  geomorphological  characteristics  studied.  All  the 
measures  under  consideration  (e.g.  Strahler’s  bifurcation,  length,  slope  and  area  ratios) 
were  within  reasonable  limits.  The  modified  model  still  presented  the  inherent  variability  ir 
network’s  growth  that  appeared  in  the  WBR  model  and  which  is  the  subject  of  study  of 
Chapter  4.  The  influence  of  diffusive  sediment  transport  processes  (e.g.  rockslide  and 
rainsplash)  was  found  to  be  more  important  in  the  modified  model  than  in  the  original  one. 
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The  influence  of  these  processes  was  shown  to  be  particularly  important  in  the  time 
evolution  of  the  hypsometric  curve.  Finally,  various  geomorphological  characteristics  (e.g. 
area-slope  renormalization,  link  length  distribution,  etc.)  were  studied  in  the  simulated 
catchments  and  their  behavior  was  found  to  be  similar  to  field  data. 

A  non-dimensional  form  of  the  modified  model,  similar  to  the  original  non- 
dimensional  WBR  model,  was  developed  and  a  new  non-dimensional  number  introduced  to 
describe  the  threshold  value  that  distinguishes  between  saturated  and  unsaturated  hillslope 
nodes.  The  purpose  of  these  non-dimensional  numbers  is  to  provide  a  way  to  meaningfully 
compare  catchments  at  different  scales. 

A  second  modification  to  the  WBR  model  was  developed  in  which  the  probabilistic 
nature  of  rainfall  was  taken  into  account.  In  the  first  modification  the  rainfall  at  the  mean 
flood  was  used  to  determine  which  nodes  were  saturated  and  which  were  not.  In  the  second 
case,  rainfalls  above  the  mean  that  saturate  and  erode  nodes  previously  considered  as 
unsaturated  were  included  in  the  model.  Differences  were  found  in  the  evolution  of  the 
hypsometric  curves.  Other  geomorphological  measures  were  not  very  different. 

An  important  conclusion  is  the  robustness  of  the  WBR  model  in  relation  to  overland 
flow  distribution.  Even  though  some  properties  of  the  behavior  change,  the  overall 
characteristics  of  the  catchment  evolutions  are  very  similar. 

The  second  part  of  this  work  studies  the  sensitivity  to  initial  conditions  in  the  WBR 
model  and  is  presented  in  Chapter  4.  Using  the  relief  difference,  a  measure  of  the  difference 
between  elevation  fields,  the  sensitivity  of  the  model  to  small  variations  in  the  initial 
elevation  field  was  studied.  This  sensitivity  is  a  manifestation  of  the  underlying  transient 
chaos  present  in  the  system.  A  variety  of  experiments,  based  on  various  tests  described  in 
the  chaos  literature,  were  used  to  analyze  the  behavior  of  the  model  under  different 
parameter  values.  Analytical  calculations  to  prove  the  existence  of  chaos  are  very  difficult 
to  develop  in  the  WBR  model  because  of  the  way  contributing  areas,  a  fundamental  term  in 
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almost  every  expression  of  the  model,  are  calculated.  Nevertheless,  the  series  of 
experiments  shown  in  Chapter  4  demonstrates  how  the  inherent  variability  of  the  resulting 
model  networks  is  related  to  extreme  sensitivity  to  initial  conditions. 

5.2  Further  Research 

Both  the  analysis  of  the  WBR  model  in  [47]  and  in  this  work  have  begun  to  look  at  a 
channel  network  evolution  model  that  keeps  the  number  of  parameters  as  small  as  possible. 
Based  on  these  analyses  a  number  of  possibilities  are  open  for  further  research.  This 
research  can  proceed  in  three  directions:  modification  and  improvement  of  the  numerical 
algorithm  that  solves  the  model;  extensions  of  the  model  to  include  different 
geomorphological  and  hydrological  hypotheses  about  the  evolution  of  the  basin;  and  field 
measurements  that  verify  and  provide  a  path  for  new  theoretical  work.  Ideally,  these  three 
branches  should  be  developed  in  parallel. 

First  of  all,  underlying  any  proposed  avenue  of  further  research  there  exists  the 
problem  of  improving  the  numerical  simulation  scheme.  The  original  code  is  time  intensive 
and  the  modified  version  is  even  more  so.  This  is  because  of  the  overland  flow  calculation 
that  has  to  take  into  account  the  saturation  deficits  at  unsaturated  nodes.  Reducing  the 
computer  time  will  allow  larger  size  domains  to  be  studied,  which  are  necessary  to  develop 
some  of  the  ideas  that  will  be  described  in  this  section. 

There  are  some  interesting  extensions  and  modifications  to  the  WBR  model.  For 
example,  the  study  of  the  effect  of  geomorphological  controls  different  from  the  tectonic 
input  used  in  the  WBR  model.  Soil  layering,  horizontal  variation  of  soil  properties  and 
limits  in  sediment  transport  produced  by  the  amount  of  material  available  for  erosion,  are 
three  possible  problems  to  study.  More  realistic  evolutions  would  probably  arise  from  these 
simulations. 
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The  tectonic  input  term  that  was  used  both  in  this  work  and  in  [47]  was  constant  in 
time  and  space.  The  effect  of  a  history  of  discrete  uplift  events  (e.g.  produced  by 
earthquakes)  instead  of  a  continuous  process,  is  also  a  possibility  to  be  studied.  Also,  the 
issue  of  erosion-deposition  and  its  relationship  with  geomorphological  configurations  like 
deltas  may  be  a  point  to  examine  carefully  with  the  model. 

The  boundary  conditions  used  in  all  the  simulations  correspond  to  idealized  walls 
around  the  domain  and  a  comer  that  was  defined  as  outlet.  Different  boundary  conditions 
would  allow  the  study  of  the  competition  for  water  and  space  between  neighboring 
watersheds  and  the  behavior  of  their  boundaries. 

Another  avenue  of  theoretical  research  comes  from  the  inclusion  of  the  probabilistic 
distribution  of  flood  rainfalls.  In  Section  3.11  the  influence  of  this  distribution  on  the 
overland  flow  production  was  studied.  Notice  that  the  indicator  function  makes  a  clear  and 
permanent  distinction  between  channel  and  hillslope  nodes.  However,  a  node  defined  as  a 
hillslope  using  the  mean  flood  rainfall  could  erode  as  a  channel  under  larger  rainfall  events. 
The  inclusion  of  this  effect  would  smooth  out  the  transition  between  channel  and  hillslope 
nodes  (Willgoose,  personal  communication).  At  this  point  sediment  transport  field  data 
linked  to  floods  as  well  as  to  the  advance  and  retreat  of  channel  heads  are  necessary. 

As  Morisawa  states  in  [26],  vegetation  plays  an  important  role  in  the  network 
expansion  process.  After  disruptive  tectonic  events  like  earthquakes,  where  vegetation  is 
sometimes  destroyed,  the  protection  of  the  soil  decreases  and  the  variation  in  elevations 
make  the  network  evolve.  As  time  passes  and  vegetation  regenerates,  the  characteristics  of 
the  process  change  and  a  temporal  variation  of  the  parameters  of  the  model  could  mimic 
this  effect. 

As  described  in  Chapter  3,  even  though  the  overland  flow  production  mechanism  is 
different  between  the  WBR  model  and  the  modified  version  developed  in  this  work,  there 
are  some  features  that  are  common  to  both  models.  Therefore,  different  simulation  tools 
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may  help  to  increase  the  simulation  speed  at  the  cost  of  losing  certain  details  of  the 
evolution  but  preserving  some  features  of  the  model’s  evolution.  Cellular  automata,  for 
example,  which  involve  discretization  not  only  in  time  but  also  in  the  variables  of  interest 
(elevation  in  this  case),  could  be  useful  in  large-scale  simulations.  These  simulations  with 
large  domain  sizes  could  be  used  in  studies  that  compare  fractal  properties  of  the  surface  or 
the  network  generated  by  the  model  with  those  measured  in  the  field.  Notice  that  the  tests  to 
measure  fractal  properties  usually  require  large  amounts  of  data  implying  very  large 
simulation  domains. 

Finally,  notice  that  the  analysis  presented  in  Chapter  3  was  based  on  certain 
assumptions  on  the  subsurface  saturation  mechanism  which  should  be  looked  at  carefully. 
Processes  like  pipeflow  and  macropore  bypassing  and  their  effect  on  the  geomorphological 
evolution  of  the  basin  may  also  be  important  in  certain  situations. 

The  third  avenue  of  research  corresponds  to  the  problem  of  field  verification  of  the 
model.  Given  that  the  timescales  of  network  evolution  are  usually  very  large,  field 
measurements  of  the  whole  process  at  large  scales  will  be  limited  except  for  a  few  special 
cases.  For  example,  settings  like  mines,  volcanic  eruptions  or  recently  deforested  zones  in 
tropical  regions  with  high  precipitation  and  easily  erodable  soils  may  be  used  for  field  sites. 
Field  experiments  should  try  to  measure  characteristics  that  enable  the  calculation  of  both 
the  nondimensional  numbers  defined  in  Section  3.8  and  the  exponents  of  sediment 
transport.  Laboratory  experiments  like  the  ones  performed  in  the  70’ s  in  the  rainfall  erosion 
experiment  facility  at  Colorado  State  University  have  been  compared  against  simulations  of 
the  WBR  model,  see  [47].  Additional  experiments  that  take  into  account  the  problem  of 
scales  in  the  diffusive  sediment  transport  processes  (like  rockslide  and  rainsplash)  can  be 
performed.  The  linkage  with  real  field  data  will  permit,  in  the  long  term,  the  study  and 
prediction  of  effects  that  change  in  land  use,  change  in  topographical  characteristics  of  a 
basin  (mines,  earthquakes,  etc.),  or  even  change  in  climate  will  have  on  the  geomorphology 
of  the  basin.  The  implications  of  Chapter  4  should  be  taken  into  account  at  this  point. 
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It  is  important  to  notice  that  even  though  the  first  studies  of  the  WBR  model  have 
been  concentrated  specifically  on  geomorphological  issues,  the  hydrological  importance  of 
the  model  should  not  be  ignored.  In  the  long  term,  more  research  could  be  directed  towards 
the  understanding  of  the  possible  links  between  the  model  parameters  and  the  hydrologic 
response  of  the  basin.  This  would  allow  conclusions  about  the  evolution  of  the  hydrological 
response  as  the  geomorphology  of  the  basins  changes  with  time. 

Finally,  in  a  more  abstract  setting,  work  with  this  model  could  also  add  some 
mathematical  insight  into  the  interesting  basin  evolution  problem.  In  general  terms  the 
problem  can  be  stated  as: 

—  =  -VnF  (5.1) ' 

dt 

where  z  is  elevation,  n  is  the  unit  vector  in  the  steepest  slope  direction,  and  F  is  the 
sediment  transport  flux  which  is  usually  assumed  to  be  a  function  of  area  and  slope.  Notice 
that  the  WBR  model  is  analogous  to  this  problem  when  some  assumptions  about  the  form 
of  F  are  added.  As  stated  in  Chapter  4,  the  mapping  that  assigns  the  contributing  area  to 
each  point  according  to  elevation  gradients  has  not  been  studied  in  detail  and  is  not  easy  to 
handle.  Up  to  now,  only  perturbation  analyses  of  the  system’s  behavior  have  been 
performed  in  [40].  A  mathematical  analysis  of  the  evolution  of  the  process,  which  is  not 
only  non-linear  but  also  spatially  variable,  is  yet  to  be  developed.  The  WBR  model  could  be 
a  useful  tool  for  this  purpose. 


Appendix  A 

Parameter  Values  of  Simulations 


SIMULATION 

TS 

TA 

TD 

TSC 

o( 

io-3 

io-8 

103 

BR643SATU10 

6.59 

17.97 

0.0 

0.92 

0.1 

SATU21 

6.59 

22.47 

0.0 

0.92 

0.1 

SATU29 

32.95 

17.97 

0.0 

0.92 

0.1 

SATU31 

6.59 

17.97 

0.0 

0.92 

0.01 

SATU32 

6.59 

17.97 

2.0 

0.92 

0.1 

SATU34 

6.59 

17.97 

0.0 

0.92 

0.1 

SATU35 

6.59 

17.97 

0.0 

0.92 

0.1 

SATU36 

6.59 

17.97 

0.0 

0.92 

0.1 

SATU37 

6.59 

17.97 

0.0 

0.92 

0.1 

SATU38 

6.59 

17.97 

0.0 

4.60 

0.1 

SATU39 

6.59 

17.97 

0.0 

0.28 

0.1 

SATU40 

6.59 

13.48 

0.0 

0.92 

0.1 

SATU42 

6.59 

7.07 

0.0 

0.92 

0.1 

SATU43 

6.59 

22.47 

0.0 

1.84 

0.1 

SATU44 

6.59 

13.48 

0.0 

0.42 

0.1 

SATU45 

6.59 

22.47 

0.0 

2.80 

0.01 

Table  A-I:  Parameters  of  Simulations 

All  the  simulations  had  m,=2,  n,=2,  m5=0.4,  ns=0.3,  TT=lxlO-6,  TC=1,  except 
SATU42  that  had  m5=0.275.  All  the  simulations  had  (36=1.0  except  SATU31  with  0.1, 
SATU37  with  10.  and  SATU45  with  0.5. 
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